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ABSTRACT 



Cebeci's viscous inviscid interaction program was applied to the analysis of 
steady two dimensional incompressible flow past a NACA 65-213 airfoil at zero angle 
of attack at a Reynolds number of 240,000. Predicted boundary layer characteristics 
were found to be quite sensitive to the choice of boundary layer transition begin and 
length. Good agreement with, the experimental results of Hoheisel ei al could be 
obtained by proper choice of both transition begin and length. 
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I. INTRODUCTION 



The prediction of the stall characteristics is of major importance in aeronautical 
engineering to determine the operating limits of an aircraft. Therefore, the 
development of reliable and accurate numerical methods for predicting separated flow 
regions is one of the most challenging problems in Computational Fluid Dynamics 
( CFD ). 

In this thesis we limit ourselves to the problem of incompressible two- 
dimensional airfoil flows. Two basic methods are available to compute viscous flows 
which include regions of flow separation. The first approach is based upon a solution 
of the full Navier-Stokes equations (or some approximate form, such as the parabolized 
Navier-Stok.es equations). This approach has the disadvantage of being very expensive 
and time consuming. The second approach is based upon the so-called Viscous- Inviscid 
Interaction Method. The outer flow is computed using the inviscid flow equations. The 
inner flow (close to the airfoil) is obtained- from a numerical solution of the Prandtls 
boundary layer equation. However, in contrast to the well-known classical boundary 
layer computations the pressure cannot be prescribed a priori, but must be found 
iteratively (ke.. by viscous-inviscid interaction). This approach has the advantage of 
being much faster and more efficient than the Navier-Stokes solutions. 

The approach chosen in this thesis is based upon the viscous-inviscid interaction 
method developed by T. Cebeci and collaborators at the Douglas Aircraft Company. In 
chapter 2 the fundamental equations are summarized . Chapter 3 is devoted to a 
discussion of the inviscid flow method, especially the so-called Panel Method first 
introduced [Ref. 4] by Hess and Smith. In chapter 4 the direct and interactive 
boundary layer methods are discussed, followed by a brief explanation of the 
turbulence model used. Finally, the computer program is explained and computed 
results are presented for the NACA 65-213 airfoil and compared with detailed 
measurements by Hoheisel et al. [Ref. 14], 



10 



II. FUNDAMENTAL EQUATIONS 



This chapter presents the equations used in our analysis, which involve the 
development of: 

1. Conservation of mass (Continuity equation! 

2. Conservation of momentum (Newton’s Second Law) 

3. Inviscid How equations 

4. Boundary layer equation 

5. Turbulence model 

With these equations we can predict the behavior of a body moving through the fluid. 
In our ca^e we are dealing with an airfoil in two dimensional, steady, inviscid and 
viscous flow. 



A. CONSERVATION OF MASS 

Let us apply the principle of conservation of mass to a small volume of space 
through which the fluid can move freely. For convenience, we shall use a cartesian 
coordinate system i.r .y.:). Furthermore, in the interest of simplicity, we shall treat a 
■2-D flow, that is. one m which there is no flow along the c -axis. Flow patterns" are the 
same for any .v-y plane. As indicated in the sketch of Figure 2.1. the component of the 
fluid velocity in the \ direction will be designated by u, and in the y direction by v. The. 
net outflow of mass through the surface surrounding the volume must be equal to the 
decrease of mass within the volume. The mass flow rate through the surface bounding 
the element is equal to the product of the density, the velocity component normal to 
the surface, and area of that surface. 

A first-order Taylor series expansion is used to evaluate the flow properties at the faces 
of the element [Ref. 3] since the properties arc a function of position. Consider the 
flow out of the volume as positive, then the net outflow of mass per-unit time is the 
summation of 
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Figure 2.1 Sketch. illustrating the velocity 
and the densitv for mass flow balance 
through a fixed volume in 2-D. 
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which must equal the rate at which the mass contained within the element decreases 



- 



Combining equation 2.1 and 2.2. dividing by A .tAy. one gets 



dp d q 

jr~ ° 



C.3) 



In vector form, equation 2.3 is 



do 

It 



T.{pV) = 0 



(2.41 



Because the pressure variations that occur in relatively low speed flow are sufficiently 
smaii. the density is essentially constant. For these incompressibe flows, the continuity 
equation becomes 



du d>' 
ds d>j 

or. in vector form V • T = 0 



B. conservation of momentum 

The equation of the conservation of linear momentum is obtained by applying 
Newton's 2nd Law [Ref. 3] where the net force acting on a fluid particle is equal to the 
time rate of change of the linear momentum of the fluid particle. As the fluid moves in 
space, its shape and volume may change, but its mass is conserved. Thus, using a 
coordinate system that is neither accelerating nor rotating, called an inertial coordinate 
system, we may write 
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F = 



m ■ 



DY_ 

Dt 



(2.7) 



The velocity V of tTfluid particle is. in general, an explicit function of time t as well as 
of its position xy. Furthermore, the position coordinates .r v v of the fluid particle are 
themselves a function of time. Since the time differentiation of equation 2.7 follows a 
given particle in its motion, the derivative is frequently termed the particle or 
substantial derivative of V, since l\xy.n and x(t),y{t). 



DV 

Dt 



where, u = 



dV dV 



= " h /,’ j- 

dx d>j 



dx_ 

dt 



and c 



dV_ 

dt 

■ - l 

dt 



Therefore, the acceleration of a fluid particle is 



DV 

Dt 



dV_ 

dt 



V ■ V)f 



(2. Si 



(2.9) 



Thus, the substantial derivative is the sum of the local, time dependent changes 
that occur at a point in the flow field and of the convection in space. When the local, 
time dependent changes are zero. oVci—0 . such flows are known as steady-state 
flows. The principal forces that act on the body are those which act directly on the 
mass of the fluid element, the body forces , and those which act on its surface, the 
pressure forces and shear forces known as surface forces. The stress system acting on an 
element of the surface is illustrated in Figure 2-2. T he stress components t acting on 
the small cube are assigned subscripts. The first subscript indicates the direction of the 
normal direction to the surface on which the stress acts and the second subscript 
indicates the direction in which the stress acts. Thus, r denotes a stress acting in the 
y direction on the surface whose normal points in the \ direction. The properties of 
most fluids have no preferred direction in space: that is. fluids arc isotropic. 

Again, the stresses on the fluid clement can be obtained from a Taylor scries 
expansion. In general, the various stresses change from point to point. Thus, they 
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Figure 2.2 Stresses acting on a 2-D element of fluid. 
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produce net forces on the fluid particle, which cause it to accelerate. The forces acting 
on each surface are obtained by taking into account the variations of stress with 
position, by using the center of the element as a reference point. 

To simplify the illustration of the force balance on the fluid particle we shall 
again consider a 2-D flow, as indicated in Figure 2-2. The resultant force in the \- 
direction. for one unit length in z is 



d d 

P J x^x A y t ( “xx ) A.'A// t ( “xv I A.r A ,y 

Ox d\j 



( 2 . 10 ) 



Where f is the body force per-unit mass in the .v direction. The most common body 
force for the How fields is that of gravity. Equation 2.10 is the left hand side of 
equation 2.". For the right hand side, combine mass term with equation 2.10 in the x 
direction 



Du ■ . 

m— = (qA.rA.y| 



l-n-.r,; 



i2.ll) 



From equation 2.10 and 2:1 1. substitute into equation 2.~ divided by A.vAv we have 



ff d 

PJx + -r df y Iy - p 



dn 

F + |r ' r) " 



(2.12) 



similarlv in the v direction 



f . d_ d 



f dr 

IJ [o7 



;r- v)r 



(2.13) 



Next, we need a relation between the stresses and the motion of the fluid. For a 
fluid at rest or for in viscid fluid motion, there is no shearing stress and the normal 
stress is in the nature of a pressure. Tor fluid particles, the stress is related to the rate 
of strain by a physical law based on the following assumptions: 
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1. Stress components may be expressed as a linear function of the components of 

the rate of the strain. The friction law for l-D flow of a Newtonian fluid is a 
special case of this linear stress/rate of strain relation, i.e., t = ( on.dy ), 

where is the fluid viscosity. 

2. The relation between the stress components and strain rate components must 
be invariant to a coordinate transformation consisting of either a rotation or a 
mirror reflection of axes, since a physical law cannot depend upon the choice of 
the coordinate system. 

3. When all velocity gradients are zero (i.e., the shear stress, vanishes), the stress 
components must reduce to the hydrostatic pressure p. 

For a fluid that satisfies these criteria, in two dimensional flow 

dn 2 

= -p+2/i— - -jt(V.F) (2.14) 

OX %j 

= =,i|r.r) |2.15| 

Oy Z 






= /' 



°JL 

[dy 



Of i ’ 



(2.16) 



with the appropriate expressions for the surface stresses, substitute into equation 2.12 
and 2.13, we have 



>[ ^ + (f • V)/i 
Ot 



f _ ££+ JL 

pfj dx dx 



„ dn 
2p — ■ 
Ox 






+ 



d_ 

dy 



, On dr , 

"< T» + T, ] 



(2.17 
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dp d 
Oy ' 0i w 



On di_ 
Oy ' dx 



+ 



o r dr 



^(v-D 



(2.1S) 
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These general differential equations for the conservation of linear momentum are 
known as the Xavier - Stokes equations. When we are dealing with incompressible and 
2-D flow, then fronTequation 2.6. V.F = 0, and the body forces are neglected, and z xx 
+ T,.,. = - 2p. Therefore equation 2.17 and 2. IS can be written as 

► V' 



du 


du 


dn 


n - — j- v 


4 - 


— 


dx 


d/J 


d l 


dv 


dv 


dr 


u hr 


u 




dx 


dlj 


dt 



1 dp 

~7 r V 

p dx 


‘ d-n 


d-n' 


. dx- ~ 


du\ 


1 dp 

— - - -L r/ 


'c> 2 r , 


d 2 v' 


t 1/ 

P d'J 


dx- ' 


dir. 



(2.19) 



( 2 . 20 ) 



These are the well known Xavier - Stokes equation for 2-D . incompressible, viscous 

flow. 



C. INVISCID FLOW EQUATION 

Inviscid flow represents an ideal How. where the effects of viscosity are zero. In 
reality, this is not true because every medium has viscosity, eventhough it may be very 
small. 

Why is the inviscid flow important in dealing with a body moving in the viscous 

fluid ? 

In 1904 L.Prandtl came up with the answer. For high Reynolds Number on a 
body moving in a viscous fluid, two regions can be distinguished. The effect of 
viscosity can be neglected outside a very thin region near the body which is called 
boundary layer. Inside the boundary layer the viscous effects are important (further 
discussion in viscous flow section). This is the reason why the inviscid flow remains 
important in the computation of fluid dynamics, eventhough it represents an ideal case. 
1. Potenial Flon 

Since the flow upstream of the body is uniform then it is also irrotational ( in 
inviscid flow). 

( = Vxf 

( 2 . 21 ) 

Consider the following vector identity , if <p is a scalar function, then 



IS 



Vx (V„*) = 0 



( 2 . 22 ) 



i.e.. the curl of the gradient of a scalar function is identically zero. Comparing equation 
2.21 and 2.22. we see that 



r = (2.23) 

Equation 2.23 states that for an irrotational flow there exists a scalar function <p such 
that the velocity is given by the gradient of (p. We denote tp as the velocity potential, tp 
is a function of the spatial coordinates, i.e.. (p = (p(.v v v). And from the definition of the 
gradient in cartesian coordinates . equation 2.23 can be written as 



Os Os 

u = - — and e = — (2.24) • 

Ox Oy 



Thus, <p has the property that its partial derivative in any direction is the 
velocity component in that direction. It follows that the existence of <p is the sole 

I • 

criterion for irrotationality. The usefulness of the velocity potential in flows of 
practical significance derives from the circumstance that, for a body in relative motion 
in an originally irrotational flow, the circulation vanishes around any contour that does 
not include the body or does not intersect the boundary layer or the wake, therefore, a 
velocity potential can be found to describe the flow everywhere outside the boundary 
layer or the wake. When a flow field is irrotational. hence allowing the velocity 
potential to be defined, there is a tremendous simplification. Instead of dealing with the 
velocity components (say. u.v and w) as the unknowns, hence requiring three equations 
for three unknowns, we can instead deal with the velocity potential as one unknown, 
therefore, requiring the solution of only one equation for the flow field and the velocity 
components can be obtained from equation 2.24. Because irrotational flows can be 
described by the velocity potential (p such flows are called potential flows. 
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2. Governing Equation for irrotational. incompressible flow : Laplace Equation. 

We have seen from equation 2.6. that the conservation of mass for an 
incompressible flowlakes the form V.V = 0. In addition, for irrotational How we have 
seen in equation 2.23 , V = V(p. Therefore, for a How that is both incompressible and 
also irrotational. equations 2.6 and 2.23 can be combined to yield. 



r- 



= o. 



or 



d 2 u d 2 v 

W 



0 



(2.25) 



Equation 2.25 is called Laplace's Equation . one of the most extensively studied 
equations in mathematical physics. 

Note that Laplace's equation is a second order linear partial differential 
equation. The fact that it is linear is particularly important, because the ’superposition 
of any particular solution of a linear differential equation is also a solution of the 

equation. For example, if <Pp (p.,. (p n represent n separate solutions of equation 

2.25, then the sunt <p = *'• <p, + + <p n is also a solution of equation 2.25. 

Since irrotational. incompressible flow is governed by Laplace's equation and 
Laplace's equation is linear, we conclude -that a complicated flow pattern for an 
irrotational. incompressible flow can be synthesized by superposition, of elementary 
flows which are also irrotational and incompressible. The singularity (or panel) 
methods presented in the next chapter are based on this idea. 

D. BOUNDARY LAYER EQUATION 

Up to this point, we have been dealing with flow outside the boundary layer, 
where viscous effects remain small. In the region within the boundary layer, velocity- 
gradients are high even with very small viscosity. Therefore, it becomes very important 
to deal with a real fluid. Before the boundary layer can be analyzed further, we have 
to know what governing equations can be used in the practical analysis. The objective 
is to predict viscous flows by means of the boundary layer method, instead of solving 
the complete Navier-Stokes equations. 

From the previous derivations, equations 2.5 , 2,19 and 2 . 20 are used here. In 
order to simplify these equations we have to make some assumptions: two- 

dimensional, steady, constant fluid properties, and no body forces. Another important 
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assumption is that the boundary layer thickness is very small compared to the length of 
the body (airfoil in this case). With these assumptions. L. Prandtl in 1904 introduced 
the "order of magnTTude" estimate into equation 2.5, 2.19 and 2.20. In order to do this, 
all linear dimensions will be referenced to the characteristic length / and all velocities 
will be referenced to U, therefore l and U can be said to have order of magnitude of 
one ( written as 0(1) ) and with the above assumptions, y will correspond to the 
boundary layer thickness (6). Then the continuity equation 2.5 can be written as 




Because 5 is very small, from the above relation we can imply that v. the velocity in 
the y direction must be very small. Therefore, the continuity equation still holds in the 
boundary layer. From the Xavier-Stokes equation in the x direction (eqn 2.19) 
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Introducing the order of magnitude in the abo\e equation, we have 
' £)•(£).« 



V 

U F 



u 



If ail above terms multiplied by 



U- 



tkeu 



P 

pU * 



Ul/u 



r-/y- 

Ul/u 



If the Revnolds Number is hish. the fourth term 



1 



= 0( zero) 



Ul/u 

l‘/6~ very large 

Ul/u large 

I- ( 

therefore — = O R. 



7-iJk 



zero = 0 ( 1 ) 



21 



So. equation 2.19 becomes: 



1 d/> d-u 

dx * d\j p dx + u dy 2 (~ 26) 



For the v direction ( equation 2.20 ) in terms of order of magnitude we obtain 
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The only term left is the pressure term, because all other terms are of higher order 

dp 



$!J 



= 0 



(2.27) 



This is very important because it tells us that the pressure across the boundary layer 
remains constant. The triplet of equations ( 2.5 , 2.26 and 2.27 ) and the boundary 
condition of zero normal and tangential velocity are known as the Prandtl's boundary 
layer equations. 

E. TURBULENT FLOW EQUATION 

Since the continuity equation (2.5) and the Navier-Stokes equation (2.19) make 
no assumptions regarding the type of flow, they are instantaneously valid in both the 
laminar and turbulent flow regimes. However, it is too difficult to deal with 
instantaneous properties in turbulent (low. Therefore we introduce the time-averaged 



properties. In our notation, prime denotes the fluctuation quantity, and bar denotes 
the mean value. 



u — n + 11 



i> = P + r 



P= P+ P 

Introducing these properties into the continuity equation and averaging, we have 



0 d 

—[a - <«') t — (r + i*') = 0 
dx dp 



Carrying out the integration term by term (details are in [Ref. 17) ). yields: 



7 — ( n — u') — -r-( r — r‘ - 
dx dp 



or 
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since — !<?) = — ( u j = - — fa) and — ( u ) — r— ( !I ) — ^ 

dx dx d x dx dx 



Therefore the continuity equation fdr turbulent flow becomes; 

-(»)- y-M = 0 

dx d;j 



(2.23) 



A similar procedure can be applied to the Navier-Stokes Equation 



d d 1 dp f d' l u 
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T-(“V) (2.29) 

dx dp 



The last two terms correspond to the normal and shear stress terms respectively, which 
we call the Reynolds Stresses or Turbulent Stresses. 



III. PANEL METHOD 



A. INTRODUCTION 

The thin airfoil theory is just what it says, it applies only to a thin airfoil at small 
angle of attack. It is not much used these days for the analysis or design of single 
element airfoils. It does give fairly good results for airfoils of 12 % thickness or less. 
On the other hand, the determination of the aerodynamic characteristics of thick, 
highly cambered, slotted surfaces, with single or multiple flaps and mutual interference 
effects among wings, fuselage, nacelles, and so forth, requires, in general, the use of 
numerical methods. 

C* 

In the 1960's Hess and Smith at McDonnell Douglas introduced a method, the 
so called Panel Method [Ref. 4j, as a numerical approach for 2-D flows which can be 
extended to 3-D potential flow problems. Such methods are called panel methods 
because the body surface is approximated by a collection of panels . There are a 
number of ways to set up the panel method. To begin with, there are choices even as 

o 

to the'type of singularity used, sources, doublets, vortices or a combination of source 
and vortex distributions. 

Panel methods as a numerical approach for predicting forces and moments acting 
on the body gave goo’d agreement with the reliable published data. The application of 
the panel method requires that the problem can be formulated such that 

1. the body can be represented by a closed polygon of a finite number of elements, 
called panels connected by nodes. 

2. the flow tangency condition is satisfied in the middle of the panels (control 
points) to avoid the inaccuracies of thin airfoil theory. 

3. the singularity distribution of each element is approximated by some kind of 
analytical function. Also, the singularities should be distributed on the body 
surface rather than on the chord line or any other line within the body. 

Before we discuss further details, it is necessary to know about the basic formulation 
for source and vortex distribution as a singularity parameter. 

1. Single Source and Source Panel 

A flowfleld where all streamlines are straight lines emanating from a source 
point, is called a source flow. On the other hand, when the flow direction is inward, it 



24 



is called sink flow. 7.1 = 0 at even. - point except at the origin where 7.1' becomes 
infinite. The source flow is irrotational at even. - point. Its velocity potential is defined 
as “ 



single source £{*•'/) — g— Inr 



( 



1) 



where A is defined as the source strength and r is the distance from the considered 
point (.r v i ) to the source. When we are dealing with non-lifting flows over a body, we 
can superimpose elementary source Hows in order to obtain a complete solution. This 
method is called source panel method. 



I 

source panel In rds ,3.2) 

. ’ o 

where * l is the panel length. 

2. Single Vortex and Vortex Panel 

The flow where all streamlines are concentric circles about a given point is 
defined as single vortex flow. The velocity along any given circular streamline is 
constant, but varies from streamline to streamline. Its velocity potential can be written 
as 



single vortex 




(3.3) 



where 9 = arctan 



'J~ Uv 



, with (j*,.. y v ) as the center of vortex 

* - J 

Let us imagine a straight line perpendicular to the page (figure 3.1). This line is a 
straight vortex filament of strength T and the flow induced in any planes perpendicular 
to the vortex of strength T. i.e., the flows in the plane 1 to the vortex filament at o 
and o' are identical to each other and are identical to the flow induced by a point 
vortex of strencth T. 



Figure 3.1 Vortex Sheet Representation. 
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Refer to Figure 3.1, imagine an infinite number of straight vortex filaments 
side by side, where the strength of each filament is infinitesimally small. Define y = 
y(s) as the strength of the vortex sheet per unit length along s, then the strength of the 
infinitesimal portion ds of the sheet is yds and the small section of the vortex sheet of 
strength yds induces an infinitesimally small velocity dv at point P(x,z) 



ds 

7 — r, _L to r 

2 “ (3.4) 



It is sometimes more convenient to use the velocity potential <p, and the increment in 
velocity potential dip induced at P[x,z) by elemental vortex yds 
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(3.5) 



For the entire vortex sheet, from a to b 







(3-6) 



Equation' 3.5 is useful in classical thin airfoil theory and equation 3.6 is important for 
the numerical vortex panel method. 



B. SOURCE AND VORTEX DISTRIBUTION 

Having introduced the basic idea of the panel method, we can use these 
singularities either alone or in combination. This method is due to Hess and Smith. 
Thus, the potential may be decomposed in a manner such that 

P = P ') + PS 4 - Pv 

(3.7) 

with <p^ being the potential of the uniform onset flow, and <p $ and <p v the potentials due 
to source and vortex distributions respectively, hence 



' /*■ _ 

T* 3 




<r(e) 



o- 



ln rds 



(3.S) 
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in which the integrations are to be performed over the body surface. Because of the 
superposition principle, this <p automatically satisfies Laplace's Equation ( see equation 
2.25) and the boundary condition at infinity. It will be the solution we seek, if cns) and 
y(s) are determined so as to meet the boundary condition of flow tangency and the 
Kutta condition (to be discussed in the next section). 

Hess and Smith assumed the vortex strength to be constant over the body 
surface and the source strength must vary over the surface. Since the Kutta condition 
involves only the trailing edge, the vortex strength can be represented by a single 
number. Thus, if one distributes on or within the body surface, vortices whose net 
strength is the correct circulation, the problem is solved if sources can be distributed 
over the body surface so as to make the total velocity field (comprised of the onset 
flow and the velocity fields due to sources and vortices) tangent to the body surface, 
regardless of how the vortices are distributed. However, the integrals in equations 3.S 
and 3.9 are hard to evaluate, even for simple forms -of the source and vortex strength, 
unless the surface on which the sources and vortices are distributed, is- a straight line. 
Thus, we select a certain number of points on the body contour, called nodes, and 
connect the nodes with straight lines, which become the panels of the method (see 
figure 3.2 ) 

We then distribute the sources and vortices on the straight line panels, so that the 
potential given by equation 3.7 can be written as: 



.v 

p = i o( i cos a + y sin o) + YL 

j = i 

In most cases, equation 3.10 still allows an exact solution of the flow problem. The 
exceptional cases are those in which the sources and vortices must be distributed 
exactly on the body surface ; to be mathematically precise, in which the potential 
cannot be continued anly tic ally across the body surface. By increasing the panel 
density, the body shape can be better approximated. This is the only major 
approximation of the panel method, one that becomes more accurate as the number of 
panels increases. 




2S 





Fig ure 3.2 Discretization: a), Actual airfoil 
b). Airlotl alter discretization. 
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Fi 2 ure 3.3 Representation ofi-th and j-th panels 
in the Panel Method. 
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To implement this method, we need a nomenclature. Let the i-th panel be 
defined as the one between the i-th and (/+ l)th nodes, and its inclination to the x axis 
be 0j. Therefore thl relation between midpoints and boundary points 



x, - 






Vi = 






and the computation of the angle 0 



B, 



— arctau 



’ 21±i 
. A', +l 




From the geometry in Figure 2.3 a relation for angle and distance between two panels 
can be obtained 



rij = \j[x> - I ])' 1 ~ r [ <J< - VjY 2 
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< j 



= arctau 



’ ( Vi - //») 
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1. Flow Tangency Condition 

The flow tangency condition in the case of no blowing or suction, is satisfied 
at the middle of each panel sometimes called the no penetration condition. To obtain 
this condition, differentiate equation 3.10 with respect to the normal vector for each 
panel. 
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AA,j ■ ■ The influence coefficient or normal velocity (normalized with F Q ) at the 
midpoint of the i-th panel due to a source distribution at the j-th panel. 

BBij • • • The influence coefficient or normal velocity (normalized with V Q ) at the 
midpoint of the i-th panel due to a vortex distribution at the j-th panel. 

, jji- . Coordinates of panel midpoint o (i-th panel 

, ,jy . Coordinates of panel midpoint of j-th panel 

(_,• length o jj-th panel 

2. Concept of the Influence Coefficient 

Equation 3.12 can be evaluated by integration. To obtain this, we ha-ve to 
relate the coordinate of a point in the j-th panel with the known coordinates of the 
boundary points and panel angles. From the geometry in Figure 3.2, one obtains 



dij_ 

dm 
d >j; 

ditj 



cos Ji 



sin Ji 



(3.13) 



(3.14) 



but, cos J3 ; = - sin 9 ; and sin p ; = cos 0 ; , therefore equation 3.13 and 3.14 become 

= - sin &i (3.15) 

= zosQi (3.1G) 

On the other hand. 

zj = Xj -r 4 j cos Oj 
!)j = Yj - 3 J sin t)j 



{ 3. IT) 
(3.13) 
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And the influence coefficients can be obtained in terms of 
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(3.19) 



AA U = - °G\ 

BB ‘‘~h '\ DF ~ CG \ ( 3.20) 



where, 



A =-(/,- .\* 7 ) cos 0j - ( ij, - !*_,) sin B } 

B = [a- Xj )' 2 +{!;,- Yj )' 2 
C = sin( Q t - B } ) 

D = cos( B, - Bj ) 

E = V 3 - A- = ( i, - ' Xj ) sin B } - ( - Y } ) cos B ] 



G — arc tan 



when i = j AA n = - and SB,, = 0 




3. Computation of Total Disturbance Velocity V 
The computation of the total disturbance velocity in x and y direction can be obtained 
by differentiating in the x and y direction. 



V = Vji + Vyj = grad f • (3.21) 
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where, 




J ~ X j ) ‘ + (i/< “ '/; )* (IS J 



- x jY 2 + {yi ~ !Jj )' 2 dsj 



. The horizontal velocity component at midpoint of i-th panel due to a unit 
source distribution at the j-th panel. 

,-L-Lfj . The vertical velocity component at midpoint of i-th panel due to a unit 
source distribution at the j-th panel. 

3Bfj . . The horizontal velocity component at midpoint of i-th panel due to a unit 
vortex distribution at the j-th panel. 

*' c 

BBf ■ • The vertical velocity component at midpoint of i-th panel due to a unit 



These influence coefficients can be obtained by integrating equation 3.12 and 
introducing equations 3.15 thru 3.18 into equation 3.12 



vortex distribution at the j-th panel. 




(3.24) 





(3.26) 
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when / = j 



.4.4f, = - £ siu Bi 

3B?;= ^cos#,- 
.4.4f, = BBf t 

BB* = -vUf, 



4. Kutta Condition 

As in all problems concerning airfoils in inviscid flows, an auxiliary condition 
needs to be invoked to ensure that a unique solution is obtained. This condition, 
known as Kutta Condition . relate to assumptions about the flow characteristics at. or at 
least in the neighbourhood of. the trailing edge. When the surface velocities are made 
equal at the midpoints of the trailing edge elements then, by the Bernoulli Equation , the 
pressures at these points are also equal, so the Kutta condition can be represented as 
the condition of zero loading in the region of the trailing edge - [Ref. 1] which is 
physically realistic. Therefore the Kutta condition can be summerized as follows: 

1. For a given angle of attack, the value of circulation T around the airfoil is such 
that the flow leaves the trailing edge smoothly. 

2. If the trailing edge angle is finite, then the trailing edge is a stagnation point. 

3. If the trailing edge angle is cusped. then the velocities leaving the top and the 
bottom surfaces at the trailing edge are finite and, equal in magnitude and 
direction. 

These are the basic principles of the Kutta condition. In the actual computation we are 
using in the code, there is no restriction whether the trailing edge angle is cusped or 
not as mentioned in item 3, but. rather, two tangential velocities in the last control 
points assume have the same magnitude but in opposite direction. 

In the numerical solutions, the actual trailing edge is not a stagnation point. 
Furthermore, it is found that the velocities at the midpoint of the trailing edge elements 
differ significantly from stagnation values: they are more likely to be closer to the free 
stream velocities. 

5. Determination of the vortex strength 

The Kutta condition is used to determine the vortex strength y. From 
equation 3-12 , set y = 1.0 and obtain the expression of <r. in terms of the influence 
coefficient 
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(3.2S) 
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This equation can~be solved by using Gaussian elimination. The result can be 
substituted into equation 3.19 and 3.20 . By letting /= 1 and /' = A’ for the trailing edge 
panel, we can compute horizontal and vertical velocities for i-th and X-th panel as ]' r j. 
]\,j. I' r Y and 1 ,.y. Introducing these values into the Kutta condition, produces a 
quadratic equation in y. This method is not the best way to satisfy the Kutta condition 
in steady How calculation, but by using this method the code can be extended to the 
unsteady case. 
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a and P, are the result obtained from equation 3.2S. From equation 3.29. there are 
two values of vortex strength but only one value is used in the further calculation. 
These values relanr to the angle of attack.. If the angle of attack is positive, use the 
negative value of y. If angle of attack is negative use the positive value. 

6. Determination of the source strength 

Once vortex strength has been determined using equation 3.12. the source 
strength <? can be obtained from 



k(= -KJ + [Jj] 



(3.30) 



7. Calculation of 'on bod} - ' velocities 

The velocity at midpoint of the i-th panel can be obtained by a spatial 
derivative of the velocity potential in the tangential direction. 



, .v S 

V, = 72 AT,j<7j + 7 72 BT,j + cos( 0i - a) 

j = i 

where AT tJ = ^\-CG- l -DF\ 



(3.31) 

(3.32) 

(3.33) 



BT.j=^-DG^\cF\ 



(3.34) 



when i = j .-IT,, = 0.0 



(3.3 5) 



3T„ = 



1 

o 



(3.36) 



AZj .is the influence coefficient of tangential direction at the midpoint of the i-th 
panel due to a unit source at j-th panel, and BT,j is the influence coefficient of 
tangential direction at the midpoint of the i-th panel due to a unit vortex at j-th panel. 
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a. Calculation of Pressure Coefficients 

Once the velocity on the body has been calculated, we can easily obtain the 
pressure coefTicientrby the Bernoulli equation 



Pi ~ PO _ , T '2 

lx-’ ~ 1 V ‘ 

•2 ^ 0 ‘ 



(3.37) 



where FI is the dimensionless velocity (normalized bv F Q ) at the midpoint of the i-th 
panel. Furthermore, as soon as the pressure coefficients are known other coefficients 
such as Lift Coefficient(CL), Drag Coefficient(CD) and Moment Coefficient about the 
leading edge (CM) can be obtained by using pressure integration along the body 
contour which is - approximated by a closed polygon. 

.v 

o, = E c„|(r , +1 - i'll 

t = i 

.V 

— H — -V+i - -V, )| 
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(3.33) 



(3.39) 



.v 



CM= E c'p.-K a * i+1 - x t )xi + (r, +l - 

i=i 



(3.40) 



CD = C x cos a + C y sin a 



(3.41) 



CL = C y cos a - C x sin a 



(3.42) 



C. LINEARLY VARYING VORTEX DISTRIBUTION 

The following method is only one variation of the use of the panel method 
[Ref. 7] ; it involves representation of the airfoil by a closed polygon of the vortex 
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panels. The vortex-panel method introduced here has the feature that the circulation 
density on each panel varies linearly from one comer to the other and is continuous 
across the comer is indicated in Figure 3.4. The airfoil and wing problems can be 
solved by means of a vortex-panel distribution alone, but calculation of fuselage and 
nacelle characteristics and their interference flows dictate the use of source and. 
possibly, doublet as well as vertex-panels. For steady flows, it has been found that the 
use of piecewise constant, discontinuous distributions of voracity can lead to 
inaccuracies such as oscillating values of the vorticity on successive panels. The use of 
a linearly varying, continuous distribution of vorticity eliminates this problem. 

In the presence of uniform flow V Q at an angle of attack a and m vortex panels, 
the velocity potential at the i-tlvc ontrol point (x y.) is defined as 



U,) = lot'.-, coset 
where ■ ^ 3 ) 

j • 



is the vortex distribution which is linear along the panel and continuous across the 
boundary points. 

The panels. N in number, are assumed planar and are named in the clockwise direction, 
starting from the trailing edge, and boundary points selected on the surface of the 
airfoil, are the intersections of continuous vortexpanels. The (.V—/) values of 7 . at 
boundary points are the unknowns to be determined numerically. The condition that 
the airfoil be a streamline is met approximately at control points. The boundary 
condition requires that the velocity in the direction outward normal vector n. be 
vanishing at the i-th control point, such that 
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Figure 3.4 Replacement of an airfoil bv vortex panels 
of linearly varying vortex strength. 



Carrying out the calculation of equation 3.44 and normalizing the vortex strengths by 2 
Jt Vq we have 
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ds j = - siul $i — a) 



(3.47) 



and from differentiation and integration we get 



.v 

E (CAT.., T ' -r CAV . V;+i | = sin(<9; - a) (3.43) 



where, y' = T--V Q is a dimensionless vortex strength. The coefficients in the 
parentheses are 



1 



CAT,, = -DF t CG - cat,, 

Cy, ;} = D+ l ;f--lAC-DE\Z- 



where 



,4 = -( X; - A* ; ) cos f)j — {iji — Yj) sin 9j 
3= (x, - A' ; )* -r ( /;, - I';)" 

C = sin (Hi - tij) 

D = cos (<9, - tij) 

E = (x, - A* ; ). sin ()j - (ili- Yj) cos 9j 



G = arc ran 

0 = ( x, — A.) cos(/9; - 2/9 ; ) - ( ji - Yj) sin( 9, - 2.9 ; ) 



/,* - 2.41, 

5 

“ Eh 
3 - .4/ , . 



The expressions in the parentheses on the left side of equation 3.48 represent the 
normal velocity at the i-ih control point induced by the linear distribution of vortices 
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on che j-ih panel, called influence coefficient of normal velocity (except for i = m-r /). 
When i = j, the coefficients have simplified values 



which describe the self-induced normal velocity at the i-th control point. 

I. Kutta Condition 

In order to apply the Kutta condition in the theoretical analysis, we need to be 
more precise about the nature of the flow at the trailing edge. It was mentioned . before 
that the trailing edge can have a finite angle or it can be cusped. The statement of the 
Kutta condition in terms of the vortex sheet is y {TE) = where V N and Vj are 

the tangential velocities at the upper and lower side of the trailing edge. However, for 
the finite-angle trailing edge, the velocities at the upper and lower surface have same 
magnitude. Thus the Kutta condition becomes y ( TE) = 0 , 



For computer programming, equation 3.48 can be arranged such 'that the- equation is 
easy to formulate in matrix form 



After combining equation 3.49 and 3.50 they are sufficient to solve for the (.V-f ■ /) 
unknown y'. values, in which the influence coefficients can be classified as 
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From the above equation, it can be expressed in matrix form 
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To obtain the unknown */'. values, the above expression can be solved by using the 
Gaussian Elimination method or any other linear equation algorithm. 

2. Calculation of the 'on body' velocities 

The unknown vortex strength having been determined, we can proceed to 
compute the velocities and pressures at every’ control point. Recall the no penetration 
condition at the control point. Hence the disturbance velocity induced at control points 
is only the tangential velocity. The velocity can be obtained by differentiating equation 
3.-- in the direction of the tangential vector on the i-ch panel, hence 



.. _ d , . 

di, 



(3.52) 



r, = coslft - o) + E ( CTujt'j + cn u l) +l | .i = 1.2.3 N (3.531 

;= i 



in which the coefficients in the parentheses are defined as 

1 



CT ltJ =-CF-DG-CT 2lJ 



1 PF 



G 



CT 2iJ = C+ -— + {AD-CZ\- 
1 



when i = j CT U , = CT lt . = -r 



The expression in the parentheses following the summation symbol has the physical 
meaning of the tangential velocity at the i-ch control point induced by the vortices 
distributed in the j-th panel. Equation 3.53 can be simplified further, to facilitate 
computer programming. 
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N 



( 3 . 54 ) 



y ti = ««(*,- - ft) -r E AZjt'j i = 1.2.3 

j=i 

aud AT, j is defined as .47,1 = C7,i 

.47, y = CT i%J + CT 2ij j = 2.3 N 

•‘U’l.V+i = CT-i iS 

in convenient matrix, the tangential velocity can be written as 
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3. Calculation of Cp.CL.CD and CM 

After the velocities at control points have been determined, one can obtain the 
pressure coefficients at the i-th control point using the Bernoulli equation 3.37 . 
Similarly, other coefficients can be obtained easily by pressure integration along the 
body contour. 



D. DISCUSSION OF THE PROGRAM PANEL 

The Panel program used in this analysis consists of two FORTRAN programs 
which involve the implementation of source and vortex distributions, and the 
implementation of linear vortex distributions. Only the first program is included in this 
thesis (see Appendix A) as a sample program. The panel methods used here are 
basically the same as the Hess and Snuth method [Ref. 4] where the source is located 
on the mid-point of each panel, constant along the panel but different from panel to 
panel . On the other hand, the vortex is considered constant for all the panels. To 
satisfy the Kutta condition at the trailing edge, the velocity components are computed 
by differentiating the velocity potential in the horizontal and vertical direction. And 
introducing these values into the Kutta condition, we can solve for the unknown vortex 
strength and subsequently for the source strength. 
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.5 Source and Vortex Panel Method. 
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In the second method, only vortices are used as singularity distributions. The only 
difference here is. that the vortices are distributed linearly along the panel and 
continuous from one panel to the other. 

1. Input Data 

The input data for the program PANEL must be arranged in the following 

order: 

1. Header card. The header card consists of the input for the number of 
coordinate points (column 1-10, integer) and input for angle of attack (column 
1 1-20 real) 

2. Coordinates of the airfoil. The arrangement for the body coordinates must be 
inputed in the following sequence: start from the trailing edge, progress on the 
lower surface to the leading edge, return through the upper surface and finish at 
the trailing edge, so that the trailing edge coordinate will be accounted twice. 

2. Program output 

The output from this program (see Appendix B) can be arranged as follows: 

TABLE 1 

DESCRIPTION OF THE PROGRAM OUTPUT 

1. PNL(I) : is the panel number 

2. X(I) : x coordinate of control points (mid-points) 

3. Y(I) : y coordinate of control points (mid-points) 

4. VEL(I) : is the dimensionless velocity for each control point 1 ' 

(total velocity divided by free stream velocity) 

5. GAMMA(I) : vortex strength) must be the same for the each panels) 

6. SIGMA(I) : source strength for each control point 

7. CP(I) : pressure coefficients for each control point 

8. CL. CD, CM : coefficient of lift, drag and moment respectively 



E. DISCUSSION OF THE RESULTS 

Figure 3.6 shows some results from the present computation for various angles of 
attack. In general, both methods give very good agreement in pressure distribution, 
although there are slight differences in the region of minimum pressure and at the 
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Figure 3 6 Comparison of pressure distributions on the N ACA _jOl . airfoil u.ing 
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trailing edge. In the leading edge neighborhood , the source strength may give a 
significant effect in the region where the velocity changes rapidly. On the other hand, 
the different implementation of the Kutta conditions at the trailing edge also give slight 
differences in the pressure distribution. Recall that the first method using both 
velocities at the first and last panels are the same in magnitude but in the opposite 
direction, while the second method is using the condition of zero vortex strength at the 
trailing edge. 

From the above conditions, the conclusions can be obtained from the present 
computations. Both methods, can be used to predict forces and moments around the 
airfoil as long as the fluid assumed remains inviscid and is not changing with time. 
Experience dictated that the second method gave less execution time than the first 
method might be true because in the second method solve only one unknown (y) 
instead of tw r o (y and <r) for the first method. On the other hand, the second method is 
more complicated in the numerical ’formulation than the first one, because of the 
varying vortex strength throughout the airfoil contour. 
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IV. VISCOUS FLOW METHOD 



In fluids, momentum is transferred by internal stresses, namely the hydrostatic 
pressure and the viscous stresses. When fluids are affected only by pressure and not by 
viscous stresses, their behavior is relatively easy to predict by standard inviscid flow 
methods as described in chapter 3. 

A % r ariation of velocity in the direction normal to the direction of the velocity 
itself is called a shear. Especially in high Reynolds Number flows this shear layer is 
very thin. 

The most common type of a shear layer is the boundary layer between a stream and a 
solid surface. On the solid surface, the fluid velocity is reduced to zero ( no slip 
condition ), but there is no direct constraint on the velocity gradient at the surface. At 
the outer edge, the velocity tends asymptotically to the free stream values. 

As mentioned before, in 1904 L. Prandtl came up with his well known theory of 
boundary layers. Prandtl's hypothesis divides the flowfield past a body into two 
separate regions, namely: 

1. The region very close to the body where viscous effects are important. 

2. The remaining region where inertia terms are more dominant than viscous 
terms, so that this region can be treated as inviscid flow. 

These assumptions allow us to deal with the parabolic boundary layer equations, 
instead of the elliptic Navier-Stokes equations. The prior experience indicates that 
parabolic equations can be solved very rapidly and efficiently. Numerically, the change 
of characteristics means a change from a field procedure to a marching method, which 
integrates the boundary layer equation for given initial conditions step by step, 
proceeding in the downstream direction. Depending on the boundary conditions, 
boundary layer methods fall into three types: 

1. The direct boundary layer method. This method employs the 'no slip 
condition',' requiring zero normal and zero tangential velocity at the surface, 
and a prescription of the external velocity at the edge of the boundary layer. 

2. The inverse boundary layer method. The prescription of wall shear or 
displacement thickness replaces the above edge boundary condition. 
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3. The interactive boundary layer method. The edge boundary condition 
prescribes a combination of displacement thickness and external velocity. 

Further discussion tvill focus on the first and the third method, since the computer 
code uses those. 



A. DIRECT BOUNDARY LAYER METHOD 



Direct methods are used in the region where the viscous effects remain small. The 
current code applies this method in the region near the leading edge. Direct methods 
allow the generation of initial conditions at the stagnation point and the efficient 
integration around the leading edge. The numerical approach features a finite 
difference method, which recasts the continuity and momentum equation as a system 
of linear algebraic equations [Ref. 13.] To begin with, we consider 2-D. steady flows of 
incompressible fluids, described in a curvilinear coordinate system with x directing 
along the airfoil surface and y perpendicular to the airfoil surface. The velocity 
components u and v shall be determined such that they satisfy anvwhere in the flow 
field the continuity equation (4.1) and the momentum equation (4.2) 



dn 

dx 

dn 




dn _ du. 
dy dx 



d 

v—(b— ) 

d !j dy 



(4.1) - 

(4.2) 



where 3 boundary conditions are needed at the boundaries of the flowfield, 

</ = 0 : u(x. 0) = 0. c( x. 0) = 0 

U = !Je • '<( x. y e ) = u,[x) 



with b * 1 + v v. These equations are referred to as boundary layer or thin shear 
layer equations. To solve these equations, it is convenient to introduce a stream 
function ( u = c\y dy and v = - <?vj/ d.x ). which reduces the number of dependent 
variables. Since the stream function automatically satisfies the continuity equation, 
only the momentum equation is left 



dc d - 1 . ’ d d‘i.' 

dy dxdy dx dy~ 
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dx Oy dy 



(4.3) 
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This equation is subjected to the Falkner-Skan transformation, which scales the normal 
coordinate y and the stream function vj/ with reference to the external velocity. 



n - \/ — y 

W U£ 



/(--• n) = 



/ u.vs 



*■( u) 



With these transformations, the momentum equation becomes 




and the boundary conditions are 

'1 = 0 : 0) =o. 0) = Q 

n = 'l s : /'(-•• r, e ) = l 



where m is defined as a dimensionless pressure gradient parameter ( m = (jr u )(du dx) 
). and prime denotes differentiation with respect to rj. This is a third order partial 
dilferential equation. The solutions of this PDE are called non-similar flows, because 
they depend on both x and \\. In contrast, if the right hand -side of equation 4.4 
vanishes, and therefore the solutions depend on rj only, they are known as similar 
Hows. 

1. The Box Method 

One of the most flexible methods in solving non-linear differential equations is 
the box method developed by Keller [Ref. 13] in 1970. The basic steps of the box 
method are the conversion of the governing equations into a first order system, the 
discretization of the differential equation by using central differences and two point 
averages, the linearization, and the solution of the resulting algebraic system. The 
introduction of two additional dependent variables U and V converts the third order 
momentum equation (4-4) into a first order system 



f = U 



(4.5) 



u' = r 



(4.6) 



51 



(bVy-r 



m — 1 



(4.7) 



■fV+m( 1 





and the boundary conditions become 

'/ = 0: U( x. QJ = 0. /( x. 0) = 0 

’1 — ’It '■ U[x. ;/,) = ! 

Instead of dealing with continuous functions f, l\ and V, we use a set of discrete 

values of these flow properties. Let the solution domain 0 <•'<*[ 0 < i) < :jj be 

covered by a rectangular mesh (Fi 2 ure 4.1) 

Xi = 0. x, — — k, with 2 < i < I 

r/x — 0. r]j = t)j-l -r hj with 2 < j < J and rfj = q e 

We approximate all quantities whether it is a function or its derivative or a parameter 

like 'b\ in terms of nodal values and coordinate of the network. The stream function 

and its first and second derivatives with respect to r| are abbreviated at the nodes of 

network by 

l/(A-.;ii). U[xi. = {/;. U). v;} 



The solution of the parabolic boundary layer equation at a certain streamline position, 
say .r.. depends solely on the solution of upstream positions say x- m j . x-_ . . . while 

no downstream influence has to be considered. The overall solution can be obtained 
step by step with the calculation propagating from the stagnation point into the 
downstream direction. The advantage of using first order equations and central 
differences is that we can reduce the domain of dependence from all upstream x- 
stations to the immediate preceding one, hence one step of the solution procedure 
writes the governing equations for a column of net rectangles (boxes) residing in the 
subdomain 

A-* < x < x, and 0 < // < 17 j 

and solves subsequently for the nodal values of the downstream face of the rectangular 
shaped subdomain. The x-station being currently solved holds therefore the superscript 
while denotes the known flow properties of the adjacent upstream location. As 
indicated by the term central differences the equations are satisfied midway between the 
nodes. 

The two ODE 's are centered about // ; _i/ 2 ) 

And the PDE is centered about ( -ft-i/e- 'Ij-i /? ) 
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Figure 4.1 



Net rectangle for finite difference approximation. 



So equations 4.5 .4.6 and 4.7 can be written in terms of finite differences 



hj 2' J 



- IT* i . 

— L = _M-‘ X ;■* ) 

h. •>' J • k ;-w 



(4.3) 

(4.0) 



(ftn — 



- 1/2 



l'i 



— / i-i r. 






77 ‘ — 77"- 1 r< /i-i 

rr‘-l/2 1/2 ’-1/2 ..,_]/•> J /._[ /2 

c ^-v2 *; ' 2—i/2 - — *; 



:4.io) 



and the boundary conditions take the form 

^ = o. /i = o 

=1 ft 

The subdomain under consideration consists of J-l net rectangles, the flow quantities 
of each being related by a momentum equation. The equations 4.S and 4.9 link the 
dependent variables to their ^-derivatives. 

2. Nekton's Method. 

U nfortunately, the unknowns appear in nonlinear combinations. Therefore we 



introduce AVu ton 's Method to solve this nonlinear system'. The solution involves an 
iterative procedure, in which the variables are linearized around their values of the 
preceding iteration 
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where 





where k denotes the iteration counter. Substitution of these values into equations 4.S, 
4.9 and 4.10, and dropping the quadratic terms in ( hUj K . i >V‘‘ K ) lead to a 
linear system in the unknowns < 6 1 y ,K 
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3 . Keller's Block Elimination Method 

Together with the boundary conditions this system is repeatedly solved until 
Sfj*. 5 Uj*. 6 Vj'* ~are small enough to be neglected. Equation 4 . 11 , 4 . 12 , and 4.13 
can be solved by Keller's Block Elimination Method. Block-tridiagonal matrices are 
composed of submatrices, called blocks, of which only those residing on main and both 
adjacent diagonals have non-zero entries. 
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The unknowns of the linear equations are the Newton iterates of the stream function ( 
fj' K ), its first derivative ( Uj^ ) and its second derivative ( Vj'^ ). This method is 
very effective and it consists primarily of two steps: 

1. The forward step eliminates the lower diagonal of submatrices. 

2. The backward step solves the remaining system from bottom to top. 



B. INTERACTIVE BOUNDARY LAYER METHOD. 

The application of the direct method is restricted to regions where the viscous 
effects remain small. Integration of the boundary layer equations will break down at 
the point of zero skin friction. To avoid this break down, we need a method that is 
able to integrate the boundary layer equations through the point of zero skin friction. 
Further, these methods are required to account for strong interaction effects, arising 
from boundary layer separation or the rapid flow acceleration downstream of the 
trailing edge, both of which cause substantial changes in the external velocity 
distribution. ’ 

In contrast with direct and inverse methods, the interactive method treats the 
external velocity and displacement thickness as unknown quantities, reflecting the 
elliptic character of the outer flows. This introduces apparently one additional 
unknown into the viscous flow problem, whose solution can be obtained by using two 
methods : 

1. The eigenvalue method, or 

2. The Mechul function method. 

The second method is being preferred here, since the first method involves non-linear 
eigenvalue problems. The edge boundary' condition of the direct problem is 
supplemented by the so-called interactive boundary condition, which relates the 
unknown external velocity with its inviscid and "displacement-perturbation-related" 
contributions. Boundary layer equations in the following constitute a system in the 
unknown functions u(x.y), v(xy) and « (. xy) 
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(4.15) 
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(4.17) 



These equations consist of continuity equation (4.15) , momentum equation in x-dircction 
(4.16) and the seemingly unnecessary momentum equation in y-direction (4.17). where 
the pressure term has been expressed in terms of the external velocity. The Mechul 
function approach assumes that the external velocity be- a function of two arguments, 
resulting in the need for the trivial y-momentum equation. The reason for considering 
uj.xy) rather than u e (x> is for purely numerical reasons, i.e., such a provision allows an 
easy setup of the finite difference equations avoiding the eigenvalue technique. 

The govering equations are complemented by proper boundary conditions. The 
velocity components u and v are required to satisfy the no-slip conditions at the surface 
of the airfoil and the horizontal component must merge smoothly into the outer flow 
at the boundary layer edge. 

U = 0 : ' u(x. 0) = 0. r(-\0)=0 

d = !J<: «(-*. -J e ) = u e (x.y e ). it e (x. y t ) = u e i\x) + - 

• ** 

with uj{x) denoting the inviscid velocity distribution and the second term, called Hilbert 
integral, approximating the perturbation velocity due to viscous effects. The interactive 
method can be applied to attached and separated How regions, while direct method 
cannot do so because of their breakdown related to the Goldstein singularity, nor 
inverse methods because of their poor convergence rates. Therefore the interactive 
methods are being preferred on the main parts of the airfoil. Only at the stagnation 
point this method cannot be applied. 

The steps which turn the partial differential equations of the interactive problem 
into a linear system of algebraic equations resemble those of the direct method, so that 
only major steps and their results will be repeated here. After the introduction of a 
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stream function (u=cy cy and v = -ciij/c.t). the equations undergo a transformation, 
which scales the normal coordinate y, stream function \\>. and the external velocity u e 
with reference to the" constant velocity u Q and the local streamwise coordinate .t 



n = 




a 



/(-.//) 

n 



—== '*•( *• u) 

f^WX 
!!>!••■ //) 

«0 



where u Q is taken as the free stream velocity. The concept of constant boundary layer 

thickness, attained bv Falkner-Skan variables with u as reference velocity, has to be 

e 

abandoned because the external velocity is unknown in the interactive calculations. 
Provided that the integration of the boundary layer equations does not start jn the 
immediate neighborhood of the stagnation point, the growth of the boundary layer 
thickness can be kept limited. In terms of these so called semi-transformed coordinates 
the boundary layer equations, written as a first order system by means of two 
additional dependent variables U and V. and the boundary conditions take the form . 
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( 4 . 21 ) 
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The discretization of the flow field follows closely the above outlined procedure of the 
direct method, covering the generation of an orthogonal grid and the introduction of 
central dilferences and two-point-averages. The overall solution proceeds in the 
downstream direction, accounting for downstream travelling disturbances only. On the 
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assumption that backflow velocities are comparatively small, a stable integration can 
be carried out by adopting the FLARE approximation ( Flugyc-Lotz and Revhner). The 
purpose of a FLARE approximation is to permit the use of a downstream-marching 
algorithm in regions of backflow. This is accomplished by setting the streamwise 
convection term u cu dx equal to zero in regions of backflow. With 



designating an "on-off switch" of the streamwise convection term, the finite difference 
equations of the interactive boundary layer problem become 



Boundary conditions are expressed in terms of nodal values, whereby the evaluation of 
the integral occuring in the interactive boundary condition involves an approximation 
in the fashion of the panel method approach, leading to 



where and c-- denote a parameter and the diagonal element of the interaction matrix, 
resulting from a discrete approximation to the Hilbert integral. Averaging as well as 






( 4 ; 25 ) 




60 



centering is supposed to obey the principle that the number of generated terms 
approaches a minimum. This entails ordinary differential equations, like the y- 
momentum equation, being centered about the middle of the downstream face and 
partial differential equations, like the x-momentum equation, being centered about the 
middle of the box. 



A balance of unknowns, which occur here as vectors in four components. 

{ f ) T * 1 confirms the principal solvability of the system. The J quadruplets of 
unknowns match with a total of 4J equations, including 2(J-I) auxiliary relations. (J-l) 
x-momentum and {J-l) y-momentum equations, each of which corresponding to one of 
the (J-l) net rectangles, and 4 boundary conditions. After linearizing this system 
around -the values of the preceding iteration ('iteration counter "k respectively 
around the solution of the adjacent upstream x-station in case of the first iteration, we 
arrive at a linear svstem in the yen-ton iterates 6 f\’ K . 
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supplemented by the two components of no-siip condition, edge and interactive 
boundary condition 
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Terms have been grouped such chat known quantities reside on the right hand sides, 
while unknown quantities appear left of the equal sign. The abbreviated coefficients in 
the momentum equation are defined by 
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and the right hand side of equation (4.2S) 
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Since the overall procedure involves a repetitive linear pattern to approach the solution 
of the nonlinear system, the linear equation solver has to be as fast as possible. Fast 
algorithms take advantage of specific matrix structures, one of which is the block- 
tndiagonal structure, which can be enforced in this method by the following 
arrangement of equations and boundary conditions : 

1. First set of equations (index j = / ) 

1. Wall boundary condition prescribing no penetration 

2. Wall boundary condition prescribing zero tangential velocity 
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3. Second auxiliary relation (linking L" and V) of the first box 

4. Trivial y-momentum equation of the first box 

2. Intermediate sets of equations (index 2 ^j^J-l) 

1. First auxiliary relation (linking/ and U) of the (j-l)-st box 

2. Momentum equation of the (j-l)-st box 

3. Second auxiliary relation (linking L" and F) of the j-th box 

4. Trivial y-momentum equation of the j-th box 

3. Last set of equations (index j= J) 

1. First auxiliary relation (linking/ and L~) of the last box 

2. Momentum equation of the last box 

3. Interactive boundary condition 

4. Edge boundary condition. 

Following these instructions equations and boundary conditions can be arranged in 
matrix-vector form. 
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( 4.30) 



where {$':*} = {Sf/. SVj'*, $W}-*} T and {r'j*} = {(rjj*. (r-)'-* . (r.j;**. 

are four dimensional vectors of the unknown Newton iterates and the 
known right hand sides, respectively. The blocks in the three diagonals of the above 
matrix ate 4;<4 and can be obtained from 
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for 2 < j < J 



for 1 < j < J - 1 



The components of the right hand side vectors can be calculated from the following 
formulae 
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momentum equatiou 



r < • — 1 



r f . K — 1 






(r 3 )-= 



• < . >C — 1 



- n 






(^r = n 7 

(/•i)l-* = 0. (r 2 )i K = 0. (r 3 )i K = 

(/■jr = 0 



for 2 < j < J 

for 2 < j < J 

for 1 < ; < J - 1 
' for 1 < ; < J — 1 
g , - c ,,-/;-*" 1 - ( i - 



The numerical solution of the above system can again be achieved by Keller's 
block elimination method, which works very much like Gauss s algorithm, but. firstly, 
matrices are eliminated instead of scalars, and. secondly, quite a few manipulations can 
be saved because of sparse occupation. 
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C. INTERACTION MODEL 



The interaction model refers to the coupling of the boundary layer and the 
external inviscid flow. In principle, the interaction consists of thickening the effective 
airfoil shape by viscous displacement, which will result in a change of the surface 
pressure. Numerically the interaction between viscous and inviscid flow regions takes 
place via the boundary conditions only, which is the specification of either an 
impermeable displacement surface or a nonzero wall transpiration at the original 
surface in case of inviscid flow, respectively, the prescription of either pressure, or 
displacement thickness, or a linear combination of both in case of viscous How. If the 
viscous effects on pressure remain small, then the interaction is called weak. However, 
situations, where viscous disturbances to the inviscid flow field are substantial, demand 
the application of strong interaction. In general , interaction models can be classified 
into four types, whereby the first three types provide loose coupling of viscous and 
inviscid regions: 

1. The direct interaction method combines a direct inviscid and direct viscous How 
solver. This classical approach achieves a solution by iterative matching, of 
boundary layer calculations with prescribed pressure and inviscid calculations 
with prescribed displacement thickness. The alternate treatment of pressure and 
displacement thickness leads to a local breakdown of the procedure when slight 
changes in pressure entail significant changes in displacement thickness (see 
Figure 4-2ai. 

2. In the inverse method the roles of the displacement thickness and pressure are 
interchanged. Hence the inviscid flow equations determine the displacement 
thickness distribution, which is imposed as boundary condition on the viscous 
flow calculation. The result of which is the pressure, which closes the cycle by 
being input to the inviscid flow computation. The hierarchical manner of 
solving the complete flow problem excludes this model just as the previous one 
from handling strongly interacting regions (see Figure 4-2b). 

3. The semi-inverse interaction method is composed of a direct inviscid and an 
inverse visous flow solver. Both parts of the scheme process the same input, i.e. 
displacement thickness, and generate the same output, i.e, pressure. After each 
cycle an updated displacement thickness is relaxed based on the deviation of the 
two pressure distributions, which should coincide upon convergence (see Figure 
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(a) DIRECT METHOD 



Airfoil 

V 

Geometry 




(b) INVERSE METHOD 



Airfoil Geometry 




| c ) SEMI-INVERSE 

figure 4.2 Direct, Inverse and Semi Inverse method. 
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4-2c). The existence of a definite hierarchy between boundary layer and the 
outer flow in the above mentioned methods gives rise to a limited rate of 
feedback between both regions. Strong interaction models incorporate the outer 
flow somehow in the boundary layer calculations, for example by the following 
interaction law: 

4. Simultaneous or strong interaction methods solve the boundary layer equations 
subject to an interaction law as outer boundary condition, which retrieves the 
elliptic character of the outer flow. No definite assignment of displacement 
thickness and pressure can be made to viscous and inviscid region. Rather, both 
quantities are treated as unknowns, related by the interaction law. The 
procedure emphasizes simultaneous solution for both displacement thickness 
and pressure (Figure 4. 2d). 



The current interaction law is formulated in terms of displacement thickness 5 and the 
external velocity :< e . the latter being related to pressure by Bernoulli's equation. The 
solution of the boundary layer equations in an iterative fashion makes use of an outer 
boundary condition, m which the totaj external velocity is written according to 

= «'/(-’•) -r u eS {x) ( J - :, l) 

where u c j(x) is the inviscid external velocity and u e ^ (,r) is the perturbation due to 
displacement effects. Using the thin airfoil concept, the perturbation velocity can be 
written as the Hilbert integral 




The contribution due to viscous effects can be evaluated by means of the blowing 
velocity concept. Lighthill proved that the effect of boundary layers on the outer flow 
can be represented by a surface distribution of sources. The source strengths must be 
determined such that the virtual displacement surface becomes a streamline (see Figure 

dt'ls) _ t'f Jr. 6 * ) (4.33) 

dx u«(.f) 
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( START ) 

I 



INPUT of 

1. Airfoil geometry 

2. Angle of elteck 

3. Reynolds number 



£ 



T V wll< 5‘ 



INVISCID FLOW SOLVER 

Laplace equation is solved by a 
conformal mapping. 

INPUT; Airfoil geometry and angle of attack:, 
displacement thickness distribution 
and blowing velocity distribution 
OUTPUT: Velocity distribution u,j on 
displacement body 



U,| 



DATA INTERFACE 

Switch from conformal mapping to b.l. grid. 
Interpolate external velocity at b.l. points. 



DATA INTERFACE 

1. RELAXATION of th« 
product term “u,J* " 

2. Compute BLOWING VELOCITY 
distribution V wt 

3. Switch from b.l. grid to 
conformal mapping grid. 
Interpolate displacement 
thickness distribution and 
blowing velocity distribution 
at boundary points 



U ol 



u,. <5* 



4 



VISCOUS FLOW SOLVER 

Boundary layer aquations era solved by a finite difference method. 



DIRECT B.L METHOD 

(Used close to lead, edge) 

INPUT: External velocity, 

initial guess of 
velocity profile 

OUTPUT; Velocity profile and 
b.l. characteristics 
like displacement 
thickness 6 * and 
skin friction 



INTERACTIVE B.L METHOD 

(Used anywhere alsa) 

Both external velocity and displacement 
thickness are treated as unknowns. 

INPUT: External velocity distribution 

and displacement thickness 
distribution of current (upstream) 
and previous (downstream) cycle, 
initial guess of velocity profile 

OUTPUT: Velocity profile including tha 
"viscoua' external velocity u. 
and b.l. characteristics like displ. 
thickness 6 * and skin friction 




OUTPUT of 

Praasure distribution 
B.l. characteristics for both 
upper and lower surfaca 
Velocity pronlea if required 



T 



( STOP ) 




Figure 4. 2d Viscous-Inviscid Interaction Method. 
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Figure 4.3 Concept of blowing velocity. 
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* 

where vi.r.o i denotes the 
sufficient to approximate 
approximation as follows: 



vertical velocity on the displacement surface. Since :t :s 
the correction term u £ ^ we can use the thin airfoil 



1. The upper - and lower surfaces of the airfoil will be considered as flat plates. 
This implies that the blowing velocity v(.t,0) equals half of the local source 
strength 

2. The displacement thickness is assumed to be so smail that the horizontal 
velocity components of the inviscid flow do not vary across the boundary layer 






= t*( -r. 0) = <5’) - 






(4.3d) 



The interaction region is limited to the finite region < x < «* t on either the upper or 
lower surface. The numerical implementation of the interaction law requires some 
discrete approximation of the above thin airfoil integral. Adopting the panel method 
approach, which concerns here a piecewise approximation of the continuous blowing 

9 * 

velocity d{ n.5')jdx to allow piecewise analytical integration, the integral can be written 
as a finite series 
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with [c,;. j denoting a matnx of interaction coefficients defining the relationship between 
the boundary layer thickness and the external velocity. Recalling that the boundary- 
layer calculation at a streamwise position involves Newton's iterations, the inviscid 
contribution can be included in the total external velocity of the previous Newton's 
iteration, leading to a generalized version of the interaction law 
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(4.36) 



Since displacement thickness does not belong to the dependent variables. the 

displacement thickness of the streamwise position currently being solved, must be 
expressed in terms of dependent variables to make allowance for its unknown status. 
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Fioure 4.4 Application of the Direct and Interactive Method. 



With ( (5 * J ' ^ being replaced by ,Jj ~ '•'j* !(u t )*•* and after separating known and 
unknown terms, the interaction law takes the form 




( 4 . 37 ) 



Implementation of this relation necessitates a known right hand side, which can be 
evaluated only in approximate fashion, because the term ( n e 5 a ) k ' K is not not known 
yet downstream of the current x-location. To ensure interaction also over this region, 
these terms arc taken from the previous iteration updated by some relaxation formula. 
With Ciu \/vxj m„ and g* = g*/n 0 denoting the dimensionless interaction 
coefficient and right hand side, respectively, the actually coded interaction law can be 
written in terms of semi-transformed coordinates 



This equation is being used as an outer boundary condition in the viscous flow 
computation, and relates the unknown external velocity with the unknown 
displacement thickness of the i-th boundary layer station, whereas the displacement 
thickness has been expressed by the strea injunction and external velocity. Because of 
the elliptic character of the outer flow, which has been incorporated in the boundary 
layer via the interaction law, the solution requires a global iteration, consisting of 
several viscous sweeps to be performed over both the streamwise upper and lower 
surface. 

Let the continuous function "external velocity times displacement thickness", 
denoted by D, be discretized at a finite number of streamwise positions, then the thin 
airfoil integral can be approximated by a finite scries of weighted "D s" at the very 
locations 
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(4.39) 
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The weights are the interaction coefficients c , with the first subscript indicating the 
streamwise position, where the correction term u e § is to be evaluated, and the second 
indicating the location, whose effect is accounted for. With a properly interpolated D- 
function, integration can be carried out analytically piece by piece. Provided the point 
under consideration does not fall within the limits of integration, D will be interpolated 
linearly. A piecewise linear function can be built up by overlapping triangular 
distributions, integration over which yields the coefficients whose k*i , kti-l , k* i— 1 
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(4.41) 



A linearly interpolated D-function would lead to singular integrals for the coefficients 
k- /, k- i-1 and k-i-r- 1 . Therefore D will be approximated by a polynomial of degree 
2 in the interval x-_j < i Splitting again into overlapping distributions, 

which this time are parabolic, and applying Cauchy's pnncipal value technique permits 
integration of singular integrands. The coefficient at the middle of the inducing source 
distribution is given by 
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The coefficient at the right of the inducing source distribution is given by 
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(4.43) 
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The coefficient at the left of the inducing source distribution is given by 

1 f 1 " 1 dD dZ 



•■•’l.l-cl — 



"^1+1 j J,_( "f| ~ £ 



(4.44) 



with 



. ; dD 
with — 
d Z 



D>~\ J-.-J,-! 



2(£- x t ] 



Di. t.[ 



-r I — 1 A — 1 i + i i 

D.^ 









* (-r.-ri — -r/)i-r--Ki — -r.— i 



hi 



1 x s — X,_i 
— 

I 1 



■ for x,_! < { < r,4.i 

for x, +l < Z < -r.+'i 

2 1 , 

l n 

< +1 * i J i +2 £ <4-1 I *T/ i-f2 



As indicated above, the overall solution is approached in an iterative process, in 
which alternately viscous and inviscid flow equations are being solved: 

1. Calculate the external velocity distribution u e j in an inviscid flow field by means 
of the conformal mapping method. To account for the airfoil-thickening due to 
viscous displacement, specify a nonzero normal velocity (blowing velocity) at 
the airfoil surface. 

2. March through the boundary layers of both streamwise upper and lower surface 
using the interaction law as outer boundary condition. 

3. Check convergence and quit if the criterion satisfied. 

4. If the convergence criterion is not met. prepare for another cycle. Update the 
product-term "external velocity times displacement thickness" on the basis of 
the deviation between inviscid and viscous external velocity distributions 
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(4.45) 
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where /. denotes the-countcr of global iterations. Further, compute the blowing velocity 
( 1 \ vlr - wall transpiration velocity ) distribution, which serves as boundary condition 
for the inviscid flow solution 
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and proceed with the first step. 



D. TURBULENCE MODELLING 



The preeence of additional unknown shear stresses in turbulent flows requires 
modelling assumptions to balance the number of unknowns and equations. Eddy 
viscosity models, one of which is used in the present method. [Ref. 1?.] relate turbulent 
shear stresses to mean flow quantities on an empirical basis. They draw their versatility 
from the,conven:ence of maintaining the same approach and numerical formulation for 
both laminar and turbulent Hows. According to this formulation for wall boundary 
layer flows, the turbulent kinematic eddy viscosity is defined by two separate formulas, 
one for the inner region being based on Van Driest’s approach, and the other for the 
outer region being based on a velocity defect approach 
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where A = 25i/ 
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Continuin' of the turbulent kinematic viscositv is established bv definine v,, as the 
distance from the wall where expressions for inner and outer region do agree. Since 
transition is not an instantaneous process, an intermediate status of flow is assumed 
between the laminar and fully turbulent regions. This region is taken into account by 
introducing an internuttency factor, which smears out the step-shaped change from 
kinematic to eddy viscosity. The formula here is suggested by Chen and Thyson. 
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7<r = 1 - exp 



(4.4S) 
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where R x[r denotes the transitional Reynolds Number ( u e x;\) tr i.e., Reynolds 
Number based on external velocity and streamwise location x [r at onset of transition. 
A value of 1200 was originally assigned to the empirical constant G., tr However, 
numerical experiments seem to indicate that low Reynolds number flows can be better 
modelled by values below 1200 (further discussion in the next section). The parameter 
a in the outer region formula is obtained from 
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where the nondimensional factor F denotes the ratio of total turbulence energy 
production to the shear-stress-related turbulence energy production, evaluated at the 
location of maximum shear stress • - • 
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The ratio of the time-averaged quantities above can be approximated by a function of 
Rt — ?tu/( — u'v‘)m<ix 
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Since this turbulence model is not validated for separated tlow, eddy viscosities in 
regions of backflow correspond to those of adjacent upstream station at the same r\- 
coordinate. Transition, if not available from other sources, currently is predicted by an 
empirical data correlation expressed in terms of the Reynolds Numbers based on 
momentum thickness and streamwise coordinate at onset of transition 
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( 4 . 51 ) 



R *tr = 1-174 ( 1 - j—) R>'*; 

r ^tr 

E. DISCUSSION OF THE COMPUTER PROGRAM 
1. Inviscid Flow Method 

The inviscid flow method used in the interaction method is based on conformal 
mapping which can be divided into the transformation of the region outside the airfoil 
to the region outside a unit circle and to the solution of the equations in the 
transformed plane. The transformation was achieved in three parts which together 
represent the major computational effort. 

In the first mapping, the airfoil is perturbed slightly to make the upper and 
lower surface trailing-edge points coincide. This is accomplished using a logarithmic 
mapping function and is necessary only in those cases in which the airfoil trailing edge 
has non-zero thickness. 

In the second mapping, the trailing-edge corner is analytically removed by 
applying the Karman-Trefftz mapping. 

In the last mapping, the resulting quasi-circular .shape is mapped to a perfect 
circle using an iterated sequence of applications of the fast Founer transform 
algorithm. The calculation of the flow in the transformed plane also makes use of 
Fourier analysis technique. 

The major computational effort required in the inviscid flow method is due to 
the, transformation. It is necessary to compute the transformation only once in the 
viscous inviscid flow interactions, so that the computational expense due to inviscid 
calculations can be held to a minimum. When the angle of attack increases, care must 
be taken since the displacement thickness can become fairly large, approaching 10% of 
the airfoil chord. In this case, use of the blowing velocity ( equation 4.33 ) on the 
airfoil surface produces a dividing streamline from the leading-edge stagnation point 
which approximates the edge of the boundary-layer displacement thickness. The 
inviscid How outside this dividing streamline is therefore the same as that past the solid 
body defined by this streamline. However, inside this dividing streamline, the inviscid 
flow is fictitious. In particular, the assumption that there is no pressure variation 
across this ficitious region becomes invalid as the magnitude of blowing velocity U 
increases. The approach adopted here is. therefore, to evaluate the velocity distribution 
directly on the displacement surface while still applying the blowing velocity on the 
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original airfoil surface. To avoid a discontinuity of the velocity between upper and 
lower surface, the Kutta condition is applied on the displacement surface. Requiring 
equal off-body pressures for the upper and lower trailing edge points, a quadratic 
equation can be solved for the unknown circulation. 

2. Interactive Viscous Flow Method 

In this method, the boundary layer is transformed into Falkner-Skan and 
subsequently semi-transformed coordinates. The numerical solution of the transformed 
equations is performed using an implicit finite-difference technique, which is first-order 
accurate in the streamwise direction and second order accurate in the normal direction. 
In principle, the computer program calculates: 

1. The inviscid pressure distribution for a specified angle of attack 

2. The boundary layer, including the velocity profile, skin friction, momentum and 
displacement thickness either: 

1. subject to a prescribed pressure distribution, or 

2. subject to an interactive edge boundary condition. In this case, the external 
velocity will be a part of the boundary layer results. 

These calculations are performed for a specified number of x-stations on the airfoil and 
in the wake, a number of sweeps is made on the airfoil in order to obtain a converged 
solution. The Cebeci-Smith two layer model is used here for computing the eddy 
viscosity and the transitional flow region is modelled using equation 4.4S. The 
method also employs a semi-empirical formula to predict the onset of transition. This 
criterion, equation 4.51. was proposed by Michel. The code offers two possibilities how 
to deal with transition: 

1. The loci of transition are calculated by the code. In this case Michel's criterion 
is employed to predict the onset of transition. When laminar separation occurs 
upstream of the calculated point of transition, then Michel's Criterion is 
disregarded and the onset of transition is redefined at the point of laminar 
separation. 

2. The begin of transition is specified by the user (fixed). The code has been 
modified to allow the onset of transition to be within or downstream of the 
laminar separation bubble. The previous version of the code always redefined 
transition at the point of laminar separation. 
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Fisure -1.3 Lift curves of the Wortmann FX 60-126 airfoil at Re = 700.000, 

luOO.OOO and 2000. WO 
(source of experimental resuits; Ref. IS). 
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F. DISCUSSION OF THE RESULTS 

1. High Reynolds Number Flows 

The computer program was First applied to high Reynolds-Number Rows about the 
Wortmann airfoil FX 60-126. Figure 4.5 shows results For Reynolds Numbers from 
700.000 up to 2000.000. It is seen that the lift predictions are in quite good agreement 
with experimental data. 

2. Low Reynolds Number Flows 

The subject of low Reynolds Number flows is important to a number of 
practical problems, including remotely piloted vehicles, wind turbines and propellers, 
sail planes, human powered vehicles, etc. Many boundary layer phenomena which have 
eluded accurate analytical prediction occur within this flow regime and are associated 
with laminar separation and subsequent transition of a laminar free shear layer. Figure 
4.6 shows the phenomenological features of the boundary layer on a low Reynolds 
number airfoil. 

When the Reynolds Number based on momentum thickness. Rq. is sufficiently low. the 
boundary layer remains laminar up to, including, the point of minimum pressure or 
maximum suction. At some location between the minimum pressure and the 
theoretical point of laminar separation the Reynolds Number of the boundary layer 
attains a critical value. This is referred to as the transition Reynolds Number based on 
momentum thickness. RQ tr 

The provision of reliable information of the flow characteristics of low 
Reynolds Number airfoils 'is hampered by the sensitivity of the flows, and particularly 
those involving separation and transition. 

In recent years, several theoretical studies have led to the development of 
methods for predicting airfoil characteristics at low Reynolds Numbers, but there are 
no universal methods which can accurately predict and account for a separation bubble 
in the design of efficient low Reynolds Number airfoils. 

The viscous,' inviscid interaction method was applied to the N'ACA 65-213 

airfoil at a Reynolds number of 240.000. It was found that the calculation would fail 

to converge if transition was predicted by Michel's criterion (equation 4.51) and if the 

empirical constant Gy tr was chosen to be 1200. Therefore it was decided to investigate 

the influence of the start of the transition and of the constant G., rr on the results bv 

r lr 

systematically varying both parameters. Table I of Appendix C shows the predicted 
lengths of the separation bubble which is obtained if transition is chosen to start at 64. 
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Figure 4.6 Phenomenological features of the boundary laver 
on the low Reynolds Number airfoil. 
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Fieure 4. 7 Experimental result for the NACA 65-213 (bv Hoheiscl et al.) 
\la_= o.l . Re-= 24o.ooo . 1 u_= 0.2"o . AOA =o.o dea 
la). Velocitv distribution and (b). Flow visualization. 
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65. 66. 67. 6S. 69. 70. 72. 74. or 76 % of chord and if Gy [r is chosen as 10. 20. 40. SO. 
or 120. It is seen that there are significant changes in the length of the separation 
bubble depending on the chosen parameter combination. This effect is displayed more 
clearly in Figure 4.11 and 4.12. 

An increase in G. /tr increases the transition length as well as the length of the 
separation bubble. 

Hoheisel et al. [Ref. 14] performed detailed laser- Doppler velocimetry 
measurements of this airfoil at a Reynolds number of 240.000. Their results are shown 
in Figures 4.7. 4.9. and 4.10 and it was attempted to choose that parameter 
combination which would give the best agreement with the experimental results. If the 
begin of transition is chosen at 74 % of chord and if G vtr =20. then the results shown 
in Figures 4.13, 4.14 and 4.15 are obtained. It is seen that the boundary layer profiles, 
displacement and momentum thickness distributions upstream of the separation bubble 
are in excellent agreement, whereas considerable deviations are found in the bubble. 
Finally, Figure 4.16 through 4.20 show the computed boundary layer velocity profiles 
for different values of Gy tr ranging from 10 to 50 
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Figure 4.10 Effect of variation ofXTRU for G yt =40 
on the velocity distribution. • 
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Figure 4.12 EfTcct of variation ofXTRL' (bcL’in of transition) for G vr =20 
on the distributions ol skin inction and lntcrmmency factor. ' 
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Figure 4.15 Comparison of 6 and 0 with experimental results. 
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Skin Friction Coefficient 
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igure 4.16 Boundarv lavcr profiles on the NACA 65-213 at Re = 240.000 

A OA — 0 deg and Gy^ r = 10. 
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Intermittency Factor, 7n 



Skin Fiiction Coefficient, 



AIRFOIL DATA 



Scales 



NACA 95-213. 

Angle Of Attaclc: 0.00 
Reynolds No.: 2400C0 
GGTH 20 




Symools: 
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T2... end of Transition 
SL.. Laminar Saporation 

ST. .. Turoulent Separation 
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Figure 4.17 Boundarv laver profile': on the NACA 65-213 at Re = 240.000 

‘ AOA = 0 deg and G. /tr = 20. 
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Intermittency Factor, 7 



Skin Friction Coefficient, 



AIRFOIL DATA 



Scales 
y/c — 



0.02 



NACA 95-213. 

Angle Of Attack: 0.00 
Reynolds No : 2400 00 
GGTR : 30 




Symbols: 



TI. .. Begin of Transition 
T2... End of Transition 
SL.. Laminar Separation 

ST. .. Turbulent Separation 
RA... Reattachment 
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Tieure 4. IS Boundary laver profiles on the NACA 65-213 at Re = 240,000 

' AOA = o deg and Gy tr = 30. 
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Intermittency Factor, 7n 
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Figure 4 19 Boundary laver profiles on the NACA 65-213 at Re = 240,000 
= ' AO A = U deg and G ytr = 40. 
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Intermittency Factor, 7i 



Skin Friction Coefficient, C,. *10 



AIRFOIL DATA 



Scales 




Symbols: T1... Begin of Transition 

T2... End of Transition 
SL.. Laminar Separation 
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Fieure 4.20 Boundarv Inver profiles on the NACA 65-213 at Re = 240.000 

' AOA = 0 deg and Gy [r = 50. 



96 



Intermittency Factor, 7t< 



V. CONCLUSION AND RECOMMENDATION 



Cebeci's viscous inviscid interaction program was applied to the analysis of 
steady two dimensional incompressible How past a NACA 65-213 airfoil at zero angle 
of attack at a Reynolds number of 240.000. Predicted boundary layer characteristics 
were found to be quite sensitive to the choice of boundary layer transition begin and 
length. Good agreement with the experimental results of Hoheisel et al could be 
obtained by proper choice of both transition begin and length. Further detailed 
measurements and calculations for other airfoils at low Reynolds number are 
recommended in order to further validate the predictive capability of the 
viscous inviscid interaction method. 



o 
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APPENDIX A 

FORTRAN PROGRAM 



CSNOEXT 



k 

k 

k 

k 

k 

k 

k 

k 

k 

k 



THIS PROGRAM CALCULATE THE PRESSURE DISTRIBUTIONS IN THE 
AIRFOIL AT ANY ANGLE OF ATTACK, BY USING PANEL METHOD IN 
2-D , INVISCID , STEADY FLOW. 

IMPLEMENTATION OF VORTEX AND SOUCRE DISTRIBUTIONS ARE USED 
ALL VELOCITIES AND LENGTH ARE NORMALIZED BY FREE STREAM 
VELOCITY AND CHORD LENGTH RESPECTIVELY. 

WRITTEN BY : CAPT.INDAF PHUTUT SUBROTO 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA. , OCTOBER 1986 



k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

c ' . 

COMMON /NODE/ N 

COMMON /AAA/ AA ( 100 , 100 ) , AAX( 100 , 100 ) , AAY( 100 , 100 ) , AAT( 100 , 100 ) 
COMMON /BEB/ BB( 100 , 100 ) , BBX( 100 , 100 ) , BBY( 100 , 100 ) , BBT(100 , 100 ) 
COMMON /STAR/ BSTAR( 100 ) , CSTAR( 100 ) 

COMMON /ALBE/ ALPA(IOO), BETA(lOO) 

COMMON /VEL/ VI ( 100 ) ,VXI ( 100 ) , VYI ( 100 ) 

COMMON /XY/ X(100) ,Y(100) ,XB(100) ,YB(100) 

COMMON /SSS/ S(100) ,SIGMA(100) ,SUMB(100) ' 

COMMON /ANGLE/ AN.ANG.TPI ,TH(100) 

COMMON /PRESS/ CP (100) 

C 

C INPUT DATA : #NODE , # AOA 

C 

READ (5,1) NN , AN 
1 FORMAT (I 10, F10. 4) 



PI = 


4 . *ATAN (1.0) 


ANG 


= AN/180. *PI 


VXA 


= -COS (ANG) 


VYA 


= -SIN (ANG) 


TP I = 


.5/PI 


DO 


10 I = 1 ,NN 



10 
5 
C 

kkkkkk 



READ (5,5) XB ( I ) , YB( I ) 
CONTINUE 
FORMAT (2F 10. 5) 



kkkkk kkkk kkkkkkkkkkkkkkkkkkkk kkkk kkkkkkkkkkkkkkkkkkkk 

k 
k 
k 
k 

kk kk kk kkkkk kkkk kkkk kk k kkkkkkkk kkkkkk kkkkkkkkkkkkkkkkkkkkkkk 

c 



COMPUTE THE ANGLES AND LENGHTS OF EACH PANELS, 
CALCULATE THE COORDINATE OF CONTROL POINTS 



N = NN- 1 
SS =0.0 
DO 



TH( I) 

S(I) 

SS 

15 CONTINUE 



(I) = .5*1 


;xb< 


;i) 


+ XB ( 


; i+i ) 


(IF = .5*1 


;ybi 


i) 


+ YB ( 


!i+i) 



= AT AN 2 ( YB( 1+1 )-YB ( I ) . XB ( 1+1 ) -XB( I ) ) 

= SQRT((XB(I+1)-XB(I))**2 + (YB ( 1+1 ) -YB (I ) ) **2 ) 
= SS + S(I) 



*********************************************************** 
* * 

* COMPUTE TIME INDEPENDENT INFLUENCE COEFFICIENTS * 
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FOR NORMAL , TANGENTIAL ,X AND Y DIRECTION 



* 
k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

c 



25 



DO 25 1= 1 ,N 
DO 25 J= 1 ,N 
IF( I .EQ. J' 
AA(I,J 
BE (I , J 
AAX ( I , J 
BBX ~ “ 
AAY 
BBY 
AAT 
BBT 

ELSE 



I, J 
( I,J 
I,J 
I,J 
,I,J 



THEN 
= 0.5 

= 0.0 

-0 . 5*SIN(TH< 
0 . 5*COS ( TH< 
0.5*COS TH< 
0 . 5*SIN(TH< 
0.0 
0.5 



A 

B 

C 

D 

£ 

F 

G 



(J) )*COS(TH( J) ) 
X(I)-XB(J) )**2 



■ jX(I)-XB 



SIN(TH(I j-THjJj j 



(Y(I)-YB(J))*SIN(TH(J)) 
( Y( I ) -YB ( J) ) x *2 



C0S(TH(I. 
(X(I)-XB(J))*SiN(TH(J)) 



= ALOG(l .0 + S(J)*(S(J)+2.*A)/B) 



. 5*C*F - D*G) 
. S^D^F + C*G) 



(Y(I)-YB(J) )*COS(TH(J) ) 



END IF 
CONTINUE 



AA( 


i, j) = 


TP I* 


BB ( 


i, j) = 


TPI* 


AAX 


?i,j) = 


TPI* 


BBX 


(i,j) = 


TP 1* 


AAY 


I,J - 


BBX( 


BBY 




-AAX( 


rttt i. 


(i'jJ = 


-BB ( I 


BBT 


(i,j) = 


AA(I 



) *F + SIN(TH(J))*Gj 



- COS (TH( J) ) *G j 



C 

C 

7rk7r7r7<7''X7rkkkkkk'Kkkkk k.k kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



* SETUP MATRICES, SET GAMA=1.0, SOLVE THE SYSTEM * 

71 BY USING GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING * 

* FOR TWO RIGHT HAND SIDES. * 

7r k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



c ‘ 

630 



63 

C 

C 



500 



C 

C 

C 

C 



DO 63 I = 1,N 

SUMB(I) =0.0 
DO 630 J =1 ,N 

SUMB(I) = SUMB(I) + BB ( I , J) 

BSTAR(I) = -SUMB(I) 

CSTAR(I) = -VXA * SIN(TH( I ) ) + VYA * COS(TH(I)) 

AA( I , N+l ) = BSTAR(I) 

AA(I,N+2) = CSTAR(I) 

CONTINUE 



CALL GAUSS (2) 

DO 500 I =1 ,N 

ALPA(I) = AA( I , N+l ) 
BETA (I ) = AA( I ,N+2) 
CONTINUE 



- MEET KUTTA CONDITIONS , SOLVE FOR VORTEX STREGNTH ( GAMA ) 
USING QUADRATIC EQUATION 



AX1 = 0.0 
AXN =0.0 
AY1 =0.0 
AYN =0.0 
BX1 = 0.0 
BXN =0.0 
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ooooo n. ooo 



510 

C 



C 



BY1 = 0.0 
BYN =0.0 

DO 510 J =1 ,N 

AX1 = AX1 + AAX ( 1 , J ) *ALPA ( J ) +BBX ( 1 , J ) 
AXN = AXN + AAX(N,j)*ALPA(j)+3BX(N,J) 
AY1 = AY1 + AAY ( 1 , j)*ALPA( J )+BBY ( 1 , j) 
AYN = AYN + AAY (N , J ) *ALPA( J )+BBY (N , J) 
BX1 = BX1 + AAX(1, J)*BETA(J) 

BXN = BXN + AAX(N, j)*BETA(j) 

BY1 = BY1 + AAY(l, J)*3ETA(J) 

BYN = BYN + AAY(N, j)*BETA(j) 

CONTINUE 



EE 

PP 



+ 

+ 



QQ 



AX1**2+AY1**2 - AXN**2-AYN**2 
AX1*BX1+AY1*BY1 - AXN* B XN - AYN* BYN + (AXN-AX1 ) *VXA 
+ (AYN-AY1 ) *VYA 

BX1**2+BY1**2-3XN**2-BYN**2 +2 . * (BXN-BX1 ) *VXA +2.* 
(BYN-BY1 )*VYA 



R = PP*PP - EE *QQ 
GAMA1 = (-PP + SORT (R))/ (EE) 
GAMA 2 = ( -PP - SQRT(R) )/(EE) 

IF (ABS (GAMA1 ) .GT. ABS (GAMA2 ) ) 
GAI-IA = GAMA 2 

ELSE 

GAMA = GAI-IA 1 



END IF 



THEN 



SOLVE THE SOURCE STRENGTH SIGMA(J) --- 

DO 550 J =1 ,M 

SIGMA (J)= GAMA*ALPA( J ) + BETA(J) 
550 ■ CONTINUE 

CALL CPRESS (GAMA) 



-- PRINT LIFT COEFFICIENT AND MOMENT COEFFICIENT ABOUT 

LEADING EDGE 



CFX =0.0 
CFY = 0.0 

CM = 0.0 • 

DO 110 I = 1 ,N 

DX = XB ( 1+1 ) - XB ( I ) 

DY = YB ( 1+1 ) - YB ( I ) 

CFX= CFX + CP(I)*DY 
CFY= CFY - CP (I *DX 
CM = CM + CP(I)*(DX*X(I)+DY*Y(I) ) 

110 CONTINUE 

CL = CFY*COS (ANG) - CFX*SIN(ANG) 

cd = cfx*cos(ang) - cfy*sin(ang) 

WRITE (8, 120) CL 
WRITE (8, 125) CD 
WRITE (3,130) CM 

120 FORMAT ( /// , 15X, 1 CL = 1 , F10 . 5 ) 

125 FORMAT (/ , 15X, 'CD =',F10.5) 

130 FORMAT ( / , 15X , ' CMLE =',F10.5) 

PRINT*, 1 "COMPUTATION COMPLETED" 1 
STOP 
END 

7T 7r7r*7t-k-k-k'k-k-k-k-k*-k*-k'k-k-k7r*-k'k-k**-k-k7<:*-k*-k***-k ********** 

x * 

* SUBROUTINE GAUSS * 

* SOLVE SIMULTANEOUS EQUATION WITH TWO RIGHT HAND SIDES * 

* BY GAUSS ELIMINATION WITH PARTIAL PIVOTING * 

* SOLUTIONS STORES IN COLUMNS NEQNS+1 AND NEQNS+2 OF * 

* MATRIX AA(NEQNS , NEQNS+2) . * 

* (AA) = COEFFICIENT OF AUGMENTED MATRIX * 



100 



no non noon noon noon no 



* NEONS = NUMBER OF EQUATIONS * 

* HRHS. = NUMBER OF RIGHT HAND SIDES * 

k k 

kkkkknkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SUBROUTINE GAUSS (NRHS) 

C 

COMMON /NODE/ NEONS 

COMMON /AAA/ AA( 100 , 100 ) , AAX( 100 , 100 ) , AAY ( 100 , 100) , AAT( 100 , 100 ) 
N? = NEONS +1 

NTOT = NEQNSrNRHS 

GAUSS REDUCTION 



DO 40 I = 2 , NEQNS 

SEARCH FOR LARGEST ENTRY IN (I-l)TH COLUMN 
ON OR 3ELOW MAIN DIAGONAL 



IM = 1-1 

I MAX = IM 

AMAX = ABS (AA(IM, IM) ) 

DO 10 J = I, NEQNS 

IF (AMAX .GE. ABS(AA(J,IM) ) ) GO TO 10 
IMAX = J 

AMAX = ABS (AA( J , IM) ) 

10 CONTINUE 

SWITCH (I-l)TH AND IMAXTH EQUATIONS 



IF (IMAX .ME. IM) GO TO 30 
DO 20 J = IM ,NTOT 

TEMP = AA(IM, J) 

AA(IM.J) '= AA( IMAX, J ) 

AA( IMAX, J) = TEMP 

20 CONTINUE 

ELIMINATE (I-l)TH UNKNOWN FROM ITH THRU 
( NEQNS )TH EQUATIONS 

20 DO 40 J = I, NEQNS 

R = AA( J , IM)/AA( IM, IM) 

DO 40 K = 'I , NTOT 
AA ( J , K ) = AA(J,K)-R*AA(IM,K) 

40 CONTINUE 

BACK SUBSTITUTION 

DO 70 K = NP,NTOT 

AA ( NEONS , K) = AA ( NEONS , K ) / AA ( NEQNS , NEQNS ) 
DO 60' L = 2, NEONS 
I = NEQNS -rl-L 

I? = I + 1 

DO 50 J = IP, NEONS 

50 AA ( I , K ) = AA(I,K)-AACl,J)*AA(J,K) 

60 AA ( I , K ) = AA(I,K)/AA(I,I) 

70 CONTINUE 
RETURN 
END 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



k k 

* SUBROUTINE CPRESS * 

* CALCULATE PRESSURE COEFFICIENTS AND TOTAL VELOCITY * 

* AT MID POINTS (CONTROL POINTS) * 

k k 



kkkkkkk^kkkkkkkkkkkkkkkkkkkkkTrkkkkkkkkkkkkkkkkkkkkkkkkkkkk^:^ 

SUBROUTINE C? RE S 3 ( GAMA ) 
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c 

COMMON /AAA/ AA( 100 , 100 ) , AAX( 100 , 100 ) , AAY( 100 , 100 ) , AAT( 100 , 100 
COMMON /BB3/ BB(100,100) , B3X( 100 , 100 ) , BBY( 100 , 100 ) , BBT ( 100 , 100 
COMMON /PRESS/ CP (100) 

COMMON /NODE/ N 

COMMON /ANGLE/ AN , ANG ,TPI ,TH( 100 ) 

COMMON /SSS/ S(100) .SIGMA(IOO) ,SUMB(100) 

COMMON /V EL/ VI ( 100 ) , VXI ( 100 ) , VYI ( 100 ) 

COMMON /XY/ X(100) ,Y(100) ,XB(100) ,YB(100) 

WRITE (8, 5) 

5 FORMAT (//,15X, 'NACA 23012' ) 

WRITE (3, 6) N 
WRITE (3, 7) AN 

6 FORMAT (15X, 'NUMBER OF PANELS =',I3) 

7 FORMAT (15X, 'ANGLE OF ATTACK =',F7.4) 

WRITE(3,20) 

DO 12 I = 1,N 
SUMM =0.0 
DO 13 J = 1 ,N 

SUMM = SUMM + AAT(I , J ) *SIGMA( J)+GAMA*BBT ( I , J) 

13 CONTINUE 

VI(I) = SUMM + COS(TH(I)-ANG) 

CP(I) = 1.0 - VI(I)**2 

■ WRITE (8,30) I,X(I) ,Y(I) ,VI(I) , GAMA , SIGMA ( I ) , CP ( I ) 

WRITE (6,31) X(I),-CP(I) 

12 CONTINUE 

20 FORMAT (// , 3X , 'PNL(I) ' ,3X, 'X(I) 1 ,7K, 'Y(I) ' ,6X, 'VEL(I) 1 ,5X, 

1 'GAMMA' ,5X, 'SIGMA(I) ' ,4X, 'CP(I) ' //) 

30 FORMAT (2X,I3,5X,F3.5,3X,F8.5,3X,F8.5,3X,F3.5,3X,F8.5,3X,F8.5) 

31 FORMAT (2F1 0.5) 

RETURN 

.END 

C 

C 
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APPENDIX B 
PROGRAM OUTPUT 



SOURCE AND VORTEX PANEL SOLUTION 
NACA 23012 

NUMBER OF PANELS = 34 
ANGLE OF ATTACK =12.0000 



PNL(I) 


X(I) 


Y(I) 


VEL(I) 


GAMMA 


SIGMA (I ) 


CP(I) 


1 


0.97500 


-0.00350 


-0.39112 


0.39035 


2.28747 


0.20590 


2 


0.92500 


-0.00965 


-0.88365 


0.39035 


2.31831 


0.21917 


3 


0.85000 


-0.01695 


-0.88620 


0.39035 


2.40999 


0.21464 


4 


0.75000 


-0.02530 


-0.88782 


0.39035 


2.41656 


0.21177 


5 


0.65000 


-0.03335 


-0.87216 


0.39035 


2.42538 


0.23933 


6 


0.55000 


-0.03920 


-0.85087 


0.39035 


2.43692 


0.27602 


7 


0.45000 


-0.04325 


-0.82394 


0.. 3903 5 


2.44450 


0.32113 


3 


0.35000 


-0.04470 


-0.77104 


0.39035- 


2.47951 


0.40549 


9 


0.27500 


-0.04370 


-0.71339 


0.39035 


2.50195 


0.49107 


10 


0.22500 


-0.04125 


-0.66069 


0.39035 


2.54546 


0.56348 


11 


0.17500 


-0.03735 


-0.57694 


0.39035 


2.60630 


0.66714 


12 


0.12500 


-0.03210 


-0.46811 


0.39035 


2.65463 


0.78087 


13 


0.03750 


-0.02765 


-0.36096 


0 . 3903.5 


2.64315 


0.36971 


14 


0.06250 


-0.02435 


-0.25392 


0.39035 


2.62380 


0.93553 


15 


0.03750 


-0.01985 


-0.00657 


0.39035 


2.65274 


0.99996 


16 


0.01375 


-0.01470 


0.41263 


0.39035 


2.58301 


0.82974 


17 


0.00625 


-0.00615 


1.41005 


0.39035 


. 2.67224 


-0.93825 


13 


0.00625 


0.01335 


2.61069 


0.39035 


-0.18981 


-5.81568 


19 


0.01375 


0.03140 


2.59588 


0.39035 


-1.31312 


-5.73860 


20 


0:03750 


0.04260 


2.32749 


0.39035 


-1.62967 


-4.41720 


21 


0.06250 


0.05355 


2.17313 


0.39035 


-1.84536 


-3.74423 


22 


0.08750 


0.06115 


2.07578 


0.39035 


-2.00811 


-3.30337 


23 


0.12500 


0.06310 


1.90635 


0.39035 


-2.23615 


-2.63608 


24 


0.17500 


0.07345 


1.74625 


0.39035 


-2.40177 


-2.04941 


25 


0.22500 


0.07550 


1.63801 


0.39035 


-2.46533 


-1.68307 


26 


0.27500 


0.07575 


1.55856 


0.39035 


’2.50044 


-1.42910 


27 


0.35000 


0.07345 


1.45669 


0.39035 


-2.55339 


-1.12194 


23 


0.45000 


0.06775 


1.35494 


0.39035 


-2.53751 


-0.83536 


29 


0.55000 


0.05940 


1.27629 


0.39035 


-2.61244 


-0.62891 


30 


0.65000 


0.04915 


1.20326 


0.39035 


-2.62994 


-0.45989 


31 


0.75000 


0.03720 


1.13842 


0.39035 


-2.65075 


-0.29600 


32 


0.85000 


0.02330 


1.07219 


0.39035 


-2.65253 


-0.14958 


33 


0.92500 


0.01300 


1.01371 


0.39035 


-2.57577 


-0.02761 


34 


0.97500 


0.00460 


0.89101 


0.39035 


-2.59479 


. 0.20611 




CL 


= 1 , 


.57532 










CD 


= - 0 , 


.63272 










CMLE 


= - 0 . 


.41146 
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APPENDIX C 
TABLE I 



TABLE 2 

EFFECT OF GGTR AND XTRL ON THE BLBBLE LENGTH 





GGTR= 10 


20 


40 


80 


120 


XTRL = .64 


.6567-. 6745 


.6567-. 6920 


.6567-. 7093 


.6567-. 7263 


.6387-7593 


XTRL = .65 


.6567-. 6920 


.6745-. 7093 


.6567-. 7263 


.6387-. 7429 


.6387-. 7752 


XTRL = .66 


.6745-7093 


.6567-. 7093 


.6567-. 7263 


.6387-. 7593 


.6204-. 7908 


XTRL = .67 


.6745-.7093 


.6745-. 7093 


.63S7-.-429 


.6204-. 7752 


.6204-. 8207 


XTRL = .68 


.6567-. 7263 


.6387-. 7593 


.6387 7 593 


.6204-. 8060 


.6020-. 8489 


XTRL = .69 


.6387-. 7429 


.6387-. 7593 


.6204-.7752 


.6020-. 8350 


.5834-. 8876 


XTRL = .70 


.6387-. 7593 


.6204-.7752 


.6204-.8060 


. .6020-. 8623 


.6204-. 8060 


XTRL =.72 


.6204-. 7908 


.6204-. 8060 


.6020-. 8350 


.5834-. 8752 


.5647-.9216 
.9415-. 9740 
(turb) 


XTRL = .74 


.6020^8060 


.6020-. 8350 


.5834-. 8623 


.5647-. 9319 


- 


XTRL =.76 


' .6020-. 8350 


.5834-.8623 


.5647-. 8752 


- 


- 



NACA 65-213 
Revnolds No. 
AO A 
XTRL 
XTRL 
GGTR 



= 240.000 
= 0.0 deg 

= Begin of iransition ( Upper Surface ) 

= Begin of transition ( Lower Surface fixed at 0.4 ) 
= Empirical Constant (Gy ) 
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APPENDIX D 
TABLE II 



TABLE 3 

EFFECT OF GGTR AND XTRL ON THE SHAPE FACTOR(H) 
AT POINT OF ZERO SKIN FRICTION 



. 


GGTR= 10 


20 


40 


80 


120 


XTRL = .64 


’ 2.830 


2.780 


2.726 


2.7423 


2.722 


XTRL' = .65 i 


2.7262 


2.745 


2.684 


2.7024 


2.7013 


1 

XTRL = .66 

! 


2.5S6 


2.624 


2.660 


2.691 

1 


2.682 


XTRL = .67 

, . ‘ 


2.604 


2.657 


2.5941 | 2.653 


2.78 


XTRL = . 68 


2.95 


2.478 


2.39 


2.6478 


•2.589 


XTRL = .69 

. . 


2.877 


2.776 


2.644 


2.5554 


2.5844 


XTRL =.70 


2.5524 


2.523 


2.534 


2.5152 


2.478 


i 

XTRL = .72 


2.54 


2.583 

i 


2.60 


2.487 


2.480 


XTRL = .74 


2.575 

1 


2.609 


2.692 


2.669 


- 


XTRL = .76 


2.624 


| 2.667 


2.82 


- 


- 



NACA 65-213 
Reynolds No. 
AOA 
H 

XTRL 

XTRL 

GGTR 



= 240.000 
= 0.(1 deg 
= 6 0 " 

= Begin of transition (Upper Surface) 

= begin of transition (Lower Surface fixed at 0.4) 
= Empirical Constant (Gy ) 
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APPENDIX E 
TABLE III 



TABLE 4 

EFFECT OF GGTR AND XTRU ON THE DRAG COEFFICIENT (CD) 





GGTR= 10 


20 


40 


80 


120 


XTRU = .64 


TE:. 01025 


.01015 


.010029 


.009876 


.009813 




WK:. 01064 


.01054 


.010423 


.010274 


.010225 


XTRU = .65 


TE:. 01022 


.01013 


.010022 


.009900 


.009785 




WK:. 01062 


.01052 


.010416 


.010301 


010193 


XTRU = .66 


TE:. 01010 


.010112 


.010018 


.009845 


.009741 




WK..01049 


.01050 


.010398 


.010248 


.010149 

i 


’ XTRU = .67 


TE:. 01016 


.01010 


.00926 


.009785 


.009816 




WK:.01055 


.01050 


.010316 


.010185 


.010236 


XTRU = .68 


TE:. 01016 


.009980 


.009944 


.009846 


.009865 


... i 


WK:. 01055 


.010383 


.010338 


.010246 


.010312 


| XTRU = .69 


TE:.01001 


.010015 


.009948 


.009886 


.010014 




WK:.01050 


.010405 


.010339 


.010316 


.010467 


XTRU = .70 


TE:.0101l 


.010015 


.009900 


.009964 


.010173 




\VK:.0l049 


.010403 


.010334 


.010382 


.010729 


XTRU = .72 


TE:. 01008 


.00997 


.009981 


.010034 


.010644 




\VK:.0i047 


.010352 


.010380 


.010481 


.011202 


XTRU = .74 


TE:.01004 


.01005 


.010017 


.010489 


_ 




WK:01042 


.010405 


.010411 


.010947 


- 


XTRU = .76 


TE:. 00999 


.010024 


.010222 


_ 






WK:. 01037 


.010434 


.010672 


- 


- 



NACA 65-213 

Revnolds No. 

AO A 

TE 

WK 

XTRU 

XTRL 

GGTR 



= 240,000 
= 0.0 deg 

= Cd at the trailing edge 
= Cd at the wake 

= Begin of transition (Upper Surface) 

= Begin of transition (Lower Surface fixed at 0.4) 
= Empirical Constant (Gy^) 
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APPENDIX F 
TABLE IV 



TABLE 5 

EFFECT OF GGTR AND XTRU ON THE LIFT COEFFICIENT (CL) 





GGTR= 10 


20 


40 


80 


120 


XTRL = .64 


.1550 


.1555 


.1559 


.1562 


.1560 


XTRU = .65 


.1557 

1 


.15602 

| 


.1565 


.1569 


.1571 


XTRU -.66 


.1571 .15630 


.1573 


• .1581 


.15*89 


XTRL -.67 : 


.1568* .15692 


.1588 


.1598 


.1599 


XTRU =. 68; .15~0 .1588 


.1594 


.1609 


.1630 


XTRU = .69 .1580 .15963 


.1609 


.1632 I .1693 


XTRU = .70 


.1596 


.16072 


.1625 


.1677 


.1800 


XTRL -.72 .1624 


.16464 


.1687 


.1815 


.2074 


XTRU = .74 


.1671 


.1708 


.1784 


.2045 


- 


XTRL = .76 

1 


.1740 


.1805 


.1910 


- 





NACA 65-213 
Revnolds No. 
AO A 
XTRU 
XTRL 
GGTR 



= 240,000 
= 0.0 deg 

= Begin of transition (Upper Surface) 

= Begin of transition (Lower Surface fixed at 0.4) 
= Empirical Constant (Gy ) 
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